Given that ∫∞−∞x2ϕ(x)dx<∞, where ϕ(x) is the standard normal probability density function, we can define the new pdf
How can I obtain a sample from f?
I know I could try brute-force inverse probability methods, but I was wondering if there is a more direct method.
Some guesswork suggest that X perhaps can be simulated by a suitable power-transformation of a Gamma random variable Y multiplied by a random sign to make the resulting density symmetric about zero. If Y has density fY(y)=λαΓ(α)yα−1e−λy,
then the density of X=YkI where P(I=−1)=P(I=1)=1/2 becomes
So for k=1/2 (a square root transformation), the gamma rate parameter λ=1/2 and the gamma shape parameter α=3/2, we obtain the desired fX.
An R implementation follows below. Note that this involves using
rgamma which uses "a modified rejection technique" (Ahrens and Dieter, 1982) so it is not clear if this is the most efficient method.
n <- 1e+4 y <- rgamma(n, shape=3/2, rate=1/2) x <- sqrt(y)*sample(c(-1, 1), n, replace=TRUE) hist(x, prob=TRUE, breaks=100) curve(x^2*dnorm(x), add=TRUE)