This tutorial should be completed with reference to these companion slides and is an optional prequel to the tutorial on parameter estimation by maximum likelihood model fitting discussed in the same slides.
Author: Sheila Kannappan (with prior major contributions from Kathleen Eckert, Rohan Isaac, and Amy Oldenberg)
Last Modified: June 2021
To complete this tutorial, download the code here and modify/run it (ideally in the Spyder app for Anaconda Python 3) with reference to the instructions below.
Computing the χ2 value is useful for determining whether a model is consistent with a data set within its errors. Most (astro)physicists define χ2 as Σi(Oi - Ei)2/σi2 where Oi, Ei, and σi are the Observed and Expected values and the errors. In other words, the numerator represents the actual residuals between the data and the model, and the denominator represents the expected residuals assuming Gaussian-distributed errors. (Note however that in online discussions, χ2 is generally defined with an Ei in the denominator, which represents the special case of Poisson-distributed data.)
To see how the χ2 value can serve as a test, consider that if the model is correct and the errors have been correctly estimated, χ2 ≈ df, where df is the number of degrees of freedom (number of data points N minus number of parameters in the model k). Therefore scientists often speak loosely and say that if the “reduced” Chi-squared defined as χ2/df is approximately equal to 1, then the fit is good. But let’s take a closer look. Open and work through the code you downloaded, uncommenting/modifying/adding code as needed and actually answering the questions embedded in the code as comments. Note, your modifications will mostly be small but make sure you understand the code provided.
(a) Using Monte Carlo methods, create 1000 fake data sets following the underlying functional form y=1/x for x = 1, 2, 3…30 with Gaussian random errors on y of amplitude 0.1. Each data set defines one value of χ2 by its residuals from the function y=1/x, and the 1000 values of χ2 from all of the data sets can be divided by df and binned into a histogram to show you the reduced χ2 distribution, which is a well-defined function analogous to a Gaussian or any other function. Note that y=1/x has no free parameters, so df is just the # of data points N=30.
(b) Now create 1000 fake data sets each with 300 values of x = 1.1, 1.2, 1.3… 30.9, using the same function y=1/x with the same errors of 0.1 on y. Overplot the new histogram of reduced χ2 values for N=300. Is a reduced χ2 of 1.3 equally good for both data sets? Google the functional form of the χ2 distribution on the web to understand why in mathematical terms.
(c) This exercise shows that just knowing that the reduced χ2 ≈ 1 does not tell you how good your model is. You must know the number of degrees of freedom (N-k). If you do, you can compute confidence levels by integrating the probability under the normalized χ2 distribution up to your measured χ2. Use np.argsort to do this approximately with the χ2 distributions from your Monte Carlo.
If you get really stuck, consult these solutions, but not 'til your puzzler is sore.