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

Given that $$∫∞−∞x2ϕ(x)dx<∞\int_{-\infty}^{\infty}x^2\phi(x)dx < \infty$$, where $$ϕ(x)\phi(x)$$ is the standard normal probability density function, we can define the new pdf

$$f(x)=x2ϕ(x)∫∞−∞t2ϕ(t)dt.f(x) = \frac{x^2\phi(x)}{\int_{-\infty}^{\infty}t^2\phi(t)dt}.$$

How can I obtain a sample from $$ff$$?

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 $$XX$$ perhaps can be simulated by a suitable power-transformation of a Gamma random variable $$YY$$ multiplied by a random sign to make the resulting density symmetric about zero. If $$YY$$ has density $$fY(y)=λαΓ(α)yα−1e−λy,f_Y(y)=\frac{\lambda^\alpha}{\Gamma(\alpha)}y^{\alpha-1}e^{-\lambda y},$$
then the density of $$X=YkIX=Y^k I$$ where $$P(I=−1)=P(I=1)=1/2P(I=-1)=P(I=1)=1/2$$ becomes
fX(x)=12fY(|x|1/k)|dydx|=12λαΓ(α)|x|(α−1)/ke−λ|x|1/k1k|x|1/k−1.\begin{align} f_X(x) &=\frac12 f_Y(|x|^{1/k})\left|\frac{dy}{dx}\right| \\&=\frac12 \frac{\lambda^\alpha}{\Gamma(\alpha)}|x|^{(\alpha-1)/k}e^{-\lambda |x|^{1/k}}\frac1k|x|^{1/k-1}. \end{align}
So for $$k=1/2k=1/2$$ (a square root transformation), the gamma rate parameter $$λ=1/2\lambda=1/2$$ and the gamma shape parameter $$α=3/2\alpha=3/2$$, we obtain the desired $$fXf_X$$.

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)