What is the most computationally efficient way of sampling $x$ from the distribution $|\psi_n(x)|^2$ for a given $n$, where $\psi_n(x)=\frac{1}{\sqrt{2^n n!}}\pi^{-1/4}\exp(-x^2/2)H_n(x)$ is a quantum harmonic oscillator eigenfunction and $H_n(x)$ is a physicist's Hermite polynomial?
The rejection method could be used, but care would be needed because $|\psi_n(x)|^2$ diverges at the classical turning points at large $n$, so that the proposal distribution would need appropriate peaks. The Metropolis algorithm should work, but may not be the most efficient solution given the need fully to decorrelate successive samples.