Sampling from $x^2\phi(x)$?

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

$$f(x) = \frac{x^2\phi(x)}{\int_{-\infty}^{\infty}t^2\phi(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.

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 $$f_Y(y)=\frac{\lambda^\alpha}{\Gamma(\alpha)}y^{\alpha-1}e^{-\lambda y},$$ then the density of $$X=Y^k I$$ where $$P(I=-1)=P(I=1)=1/2$$ becomes \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/2$$ (a square root transformation), the gamma rate parameter $$\lambda=1/2$$ and the gamma shape parameter $$\alpha=3/2$$, we obtain the desired $$f_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)


• This is a nice approach. The answer could be improved by also explaining the math behind it (i.e. how the OP's distribution relates to the gamma distribution, and how this leads to the sampling method you described). Apr 7, 2021 at 12:42
• @user20160 Alright, I added some explanation on how I came up with this Apr 7, 2021 at 13:11

One could consider the following alternatives to Jarle Tufto's most efficient solution:

1. namely, transform back into $$X$$ a Gamma $$\mathcal G(3/2,1/2)$$ [or equivalently a $$\chi^2_3$$] variate
rjt=function(n)sqrt(rgamma(n,3/2)*2)*sample(c(-1,1),n,rep=TRUE)

1. use a numerical inverse of the cdf$$F(x)=\int_{-\infty}^x y^2\varphi(y)\,\text dy=\int_{-\infty}^x y\,\varphi^\prime(y)\,\text dy=\Phi(x)-x\varphi(x)$$
pf=function(x)pnorm(x)-x*dnorm(x)
df=function(x)x^2*dnorm(x)#the denominator is the Normal variance, 1
qf=function(u)uniroot(function(x)pf(x)-u,int=c(-10,10))\$root
rf=function(n)apply(as.matrix(runif(n)),1,qf)

1. apply an accept-reject solution based on a Normal density with a larger variance$$^1$$ by noticing that$$f(x)=\frac{1}{\sqrt{2\pi}}x^2e^{-x^2/4}e^{-x^2/4}\le\frac{4}{e}\frac{\sqrt{2}}{\sqrt{4\pi}}e^{-x^2/4}=M \varphi(x;0,2)$$
arf=function(n,M=4*sqrt(2)/exp(1)){
z=rnorm(M*n,0,sd=sqrt(2))
z[runif(M*n)<df(z)/dnorm(z,sd=sqrt(2))/M]}


The correct fit of all three generators can be tested by the following qq-plot command:

 plot(1:1e4,pf(sort(rf(1e4))))


This range of solutions produces the following relative and respective running times

           test replications elapsed relative user.self
3 accept-reject           100    0.294    2.211    0.294
2       inverse           100   74.707  570.282   74.681
1     transform           100    0.131    1.000    0.131


which shows the superiority of the original transform$$^2$$ approach.

$$^1$$As an addendum, the accept-reject solution can be improved by optimising the break-up$$\exp\{-x^2/2\}=\exp\{-\alpha x^2/2\}\exp\{-(1-\alpha)x^2/2\}$$ as a function of $$\alpha$$, since it is easy to show that $$\alpha=2/3$$ minimises the upper bound $$M=2/\alpha\sqrt{1-\alpha} e$$. However, this finer tuning of accept-reject only brings a gain of less than 10%:

3 accept-reject   0.273    2.133     0.269
2       inverse  73.895  577.305    73.878
1     transform   0.128    1.000     0.128


$$^2$$Devroye (1986) mentions this distribution a few times in his Bible of simulation methods (p.119, p.176) as the Maxwell distribution. In particular, if $$X$$ follows this distribution and if $$U$$ is Uniform on $$(0,1)$$, $$Y=UX$$ follows a Normal distribution. Sadly, the reciprocal does not work: simulating a Normal variate and dividing by an independent Uniform variate does not return a Maxwell variate! However, inverting the joint distribution of $$(U,X)$$ into a joint distribution of $$(X,Y)$$ leads to the (conditional on $$Y$$) representation $$X=\frac{Y}{|Y|}\sqrt{Y^2-2\log(V)}\qquad Y\sim\mathcal N(0,1)\quad V\sim\mathcal U(0,1)$$ ie

rmax=function(n)sqrt((y<-rnorm(n))**2-2*log(runif(n)))*y/abs(y)


which proves to be 20% faster than the Gamma transform in rjt :

     test replications elapsed relative user.self
2 devroye          100   0.991    1.000     0.992
1   tufto          100   1.246    1.237     1.245


despite $$Y^2-2\log(V)$$ being just another representation of a Gamma $$\mathcal G(3/2,1/2)$$ variate.

First of all, it is worth noting that the scaling constant in this case is the second raw moment of the standard normal distribution, which is:

$$\int \limits_{-\infty}^\infty x^2 \phi(x) \ dx = 1.$$

$$f(x) = \frac{x^2}{\sqrt{2 \pi}} \cdot \exp \Big( -\frac{x^2}{2} \Big) \quad \quad \quad \text{for all } x \in \mathbb{R}.$$

You can sample from this density without rejection using the transformation (see below for proof):

$$X = \text{SGN} \cdot \chi \quad \quad \quad \text{SGN} \sim 1 - 2 \cdot \text{Bern}(\tfrac{1}{2}) \quad \quad \quad \chi \sim \text{Chi}(\text{df} = 3).$$

We can easily implement this transformation method in R to produce the following simulation function (which is vectorised to allow you to produce any number of simulations).

rtransnormdist <- function(n) {
CHI <- sqrt(rchisq(n, df = 3))
SGN <- sample(c(-1, 1), size = n, replace = TRUE)
SGN*CHI }



We can confirm that this produces the required density as follows:

set.seed(1)
SIMS <- rtransnormdist(10^6)
plot(density(SIMS), lty = 2, lwd = 2, main = 'Simulated Density')
curve(x^2*dnorm(x), col = 'red', lty = 3, lwd = 2, add = TRUE)


Proof of density transformation: Using the stated transformation and applying the rules for density transformations we obtain:

\begin{align} f_{|X|}(x) = f_\chi(x) \cdot \Bigg| \frac{dx}{d \chi} \Bigg| &= \text{Chi}(x|3) \times 1 \\[6pt] &= \frac{x^2 \sqrt{2}}{\sqrt{\pi}} \cdot \exp \Big( -\frac{x^2}{2} \Big) \cdot \mathbb{I}(x \geqslant 0), \\[6pt] \end{align}

which then gives the density:

\begin{align} f_{X}(x) = \frac{1}{2} \cdot f_{|X|}(|x|) &= \frac{x^2}{\sqrt{2 \pi}} \cdot \exp \Big( -\frac{x^2}{2} \Big). \\[6pt] \end{align}

This confirms the desired density function.