Show the pdf of a chi-squared distribution using python

Question

I'm trying to reconstruct the pdf of the chi-squared distribution with 3 degrees of freedom from a simulated sample. Here is my python code:

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

norm = stats.norm(0, 1)

x1 = [x * x for x in np.random.randn(1000)]
x2 = [x * x for x in np.random.randn(1000)]
x3 = [x * x for x in np.random.randn(1000)]

f = x1 + x2 + x3

plt.hist(f, 100)
plt.show()

The result I got was this.

Of course this is wrong. As shown in Wikipedia, the pdf of the chi-squared distribution with 3 degrees of freedom should go upwards first from zero and go downwards later, not something keep climbing like mine. Is there anything wrong with my code? The formula I used was as follows:

Q = x1^2 + x2^2 + x3^2

where x1, x2 and x3 are independent, standard normal random variables.

Answer 1

The '+' operator works differently on Python lists than on Numpy arrays.

f = x1 + x2 + x3

concatenates three lists into one. However, you want to add the content of the three lists element-wise, which could be done like this:

f = np.array(x1) + np.array(x2) + np.array(x3)

Answer 2

Although I tried your code and got the same result as you, if you use your 'norm' variable to generate the random values it seems to work.

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

norm = stats.norm(0, 1)

x1 = norm.rvs(size=100000)**2
x2 = norm.rvs(size=100000)**2
x3 = norm.rvs(size=100000)**2

f = x1 + x2 + x3

plt.hist(f, 60, normed=True)

# Plot the theoretical density of f
x = np.arange(0, 30, .05)
plt.plot(x, stats.chi2.pdf(x, df=3), color='r', lw=2)
plt.show()

The result I got was