Sampling from x2ϕ(x)x^2\phi(x)?

Given that x2ϕ(x)dx<, where ϕ(x) is the standard normal probability density function, we can define the new pdf

f(x)=x2ϕ(x)t2ϕ(t)dt.

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.

Answer

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
fX(x)=12fY(|x|1/k)|dydx|=12λαΓ(α)|x|(α1)/keλ|x|1/k1k|x|1/k1.
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)

enter image description here

Attribution
Source : Link , Question Author : motto , Answer Author : Jarle Tufto

Leave a Comment