Deriving $\chi^2$ density from the standard normal $Z$ density

Question

For the sake of learning, I am trying to derive a $\chi^2$ density from the standard normal $Z$ density.

The $\chi^2_k$ is a distribution of the sum of squares of $k$ standard normal variables, where $k$ is the the number of degrees of freedom.

Taking $k=1$ for simplicity (i.e. just one standard normal variable), I have tried the following in R:

library(data.table)

DT <- data.table(x = seq(0, 3.5, by=0.01)) #In a range from 0 to 3.5

Standard normal density multiplied by 2 to calculate the density of absolute $Z$ values (i.e. negative + positive); this is equivalent to a $\chi$ density.

DT[, Z_abs := 2*dnorm(x)] 

Calculate the squares:

DT[, X_squared := x^2]

And plot the density against the squares:

library(ggplot2)

ggplot(DT, aes(X_squared, Z_abs)) + geom_line() +
  scale_x_continuous(limits=c(0,3.5), breaks=seq(0,3.5,0.5))

However, this is apparently wrong as it differs from the correct $\chi^2$ density given by the built-in dchisq function (shown in red).

DT[, Chisq := dchisq(X_squared, df=1)]

ggplot(DT, aes(X_squared, Z_abs)) + geom_line() +
  geom_line(aes(y=Chisq), color="red") +
  scale_x_continuous(limits=c(0,3.5), breaks=seq(0,3.5,0.5)) +
  scale_y_continuous(limits=c(0,0.8)) +
  ylab("Density")

Where did I get it wrong ?

Answer 1

Let us first find the density of the absolute value of $X\sim N(0,1)$, $Y=|X|$, $g(y)$. For $y>0$, \begin{eqnarray*} G(y)\equiv\Pr(Y\leqslant y)&=&\Pr(|X|\leqslant y)\\ &=&\Pr(-y\leqslant X\leqslant y)\\ &=&F(y)-F(-y) \end{eqnarray*} and hence $$g(y)=G'(y)=f(y)+f(-y)$$ (note the second inner derivative equals $-1$). Since $y=|x|$ cannot be negative, we have $g(y)=0$ for $y<0$. For densities symmetric about zero (such as the normal density) we have $$f(y)=f(-y)$$ and therefore $$ g(y)=\begin{cases}0&y<0\\ 2f(y)&y\geqslant0\end{cases} $$ For the second step (the change of variable technique that you were referred to in the comments), we have $y=w(z)=\sqrt{z}$ and hence $w'(z)=\frac{1}{2}z^{-\frac{1}{2}}$. Thus, \begin{eqnarray*} h(z)&=&\frac{2}{\sqrt{2\pi}}\exp\biggl\{-\frac{(\sqrt{z})^2}{2}\biggr\}\left|\frac{1}{2}z^{-\frac{1}{2}}\right|\\ &=&\frac{1}{\sqrt{2\pi}}z^{-\frac{1}{2}}e^{-\frac{z}{2}} \end{eqnarray*} for $z>0$ and 0 elsewhere. This is the density of the chi-square distribution with $\nu=1$.

Answer 2

As pointed out in the comments, your error here is that your density transformation does not take account of the nonlinearity of the transformation. I will show you a better empirical demonstration of the distributional equivalence, where we don't attempt the transform at all, but simply compare the kernel density of simulated values of the sum with the postulated chi-squared density. I will also show you a simple proof of the distributional equivalence at issue.

---

Empirical simulation: Rather than attempting the density transformation, let's proceed by simulation, by simulating $N=10^5$ sets of $n=4$ standard normal random variables.

#Generate matrix of simulations
set.seed(1)
N <- 10^5
n <- 4
SIMS <- matrix(rnorm(N*n), nrow = N, ncol = n)

#Compute statistic of interest and its kernel density
STAT <- rep(0, N)
for (i in 1:N) { STAT[i] <- sum(SIMS[i,]^2) }
DENS <- density(STAT)

#Set chi-squared density function
CCC  <- function(x) { dchisq(x, df = n) }

#Plot the kernel density against the postulated chi-squared density
plot(DENS, xlim = c(0, 20), lwd = 2, main = 'Simulation of Density')
plot(CCC,  xlim = c(0, 20), lwd = 2, lty = 2, col = 'red', add = TRUE)

As you can see from the plot, the simulated values closely follow the postulated chi-squared density. You can easily repeat this simulation analysis for different values of $n$ if you would like to demonstrate the distributional equivalence for other values.

---

Proving equivalence in distribution: For completeness, I supply you here with a proof of the distributional result you are trying to show by simulation. The simplest way to prove this result is via moment generating functions (or characteristic functions). Let $Z_1, ..., Z_n \sim \text{IID N}(0,1)$ be a set of IID standard normal random variables and let $G = \sum_{i=1}^n Z_i^2$. Using the law of the unconscious statistician and the substitution $y = z \sqrt{1/2-t}$, for all $t < \tfrac{1}{2}$ we have:

$$\begin{align} \mathbb{E}(\exp(t Z_i^2)) &= \int \limits_{-\infty}^\infty \exp(t z^2) \cdot \frac{1}{\sqrt{2 \pi}} \exp \bigg( -\frac{1}{2} z^2 \bigg) \ dz \\[6pt] &= \frac{1}{\sqrt{2 \pi}} \int \limits_{-\infty}^\infty \exp \bigg( - \Big( \frac{1}{2} - t \Big) z^2 \bigg) \ dz \\[6pt] &= \frac{1}{\sqrt{(1-2t) \pi}} \int \limits_{-\infty}^\infty \exp ( -y^2 ) \ dy \\[6pt] &= \frac{1}{\sqrt{(1-2t)}}, \\[6pt] \end{align}$$

Thus, for all $t < \tfrac{1}{2}$ the moment generating function for $G$ is:

$$\begin{align} m_G(t) &\equiv \mathbb{E}(\exp(tG)) \\[10pt] &= \prod_{i=1}^n \mathbb{E}(\exp(t Z_i^2)) \\[6pt] &= \prod_{i=1}^n \frac{1}{\sqrt{(1-2t)}} \\[10pt] &= (1-2t)^{-n/2}. \\[6pt] \end{align}$$

This is the moment generating function for the chi-squared distribution, which demonstrates that $G$ is a chi-squared random variable. (See here for proof that the moment generating function determines the distribution.)