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.
 A: 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)
curve(x^2*dnorm(x), add=TRUE)


A: 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.$$
Consequently, your density function is:
$$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.
