8. Workshop:  Chi-Squared Test for Goodness of Data

The Chi-Squared (or (χ^2)) Test  for data is a measure of how good your experimental data (e.g. your background radiation data) is represented by a model distribution (e.g. the Poisson distribution).

The (x^2) Test

Suppose you have observational data (Subscript[O, k]) where k indicates bin number. In the case of the background radiation, k is just the number of counts (n=0, 1, 2, ...) on a given observation and (Subscript[O, k]) is the number of occurances of k counts in your data set.

ASIDE:  Sometimes you must combine your data into "bins".  For example, if you have a "hot" beta source which produces a large number of counts in say 10 seconds then it makes sense of combine some of your data points (e.g. those occurance where the number of counts is between 1000 and 1250 counts in 10 seconds).  If you did not do this your data set would contain a large number of instances of either zero or one occurances.

Corresponding to the observational data (Subscript[O, k]) you probably have some theoretical model you expect the observational data will fit.  The expected occurances (for the kth bin) is called (Subscript[E, k]).  (e.g. For the Poisson distribution you first take your experimental data of Num=50 total observations.  Compute the average which you call λ.  The corresponding Poisson distribution  P[n]=(((e^(-λ)) (λ^n))/(n!))gives (Subscript[E, n])= Num×P[n]  with n=k.)

ASIDE:   For the Poisson distribution you first take your experimental data of Num=50 total observations.  Compute the average which you call λ.  The λ you calculated from the experimental data defines the corresponding (theoretical) Poisson distribution P[k]=(((e^(-λ)) (λ^k))/(k!))gives 
(Subscript[E, k])= Num×P[k] .

Obviously you do not expect perfect agreement (because of the statistical uncertainty in the experiment) between you observational data (Subscript[O, k]) and you theoretical model (Subscript[E,  k]).  ((Subscript[O, k])-(Subscript[E, k])) is a measure of the of the deviation of the observation of the kth bin from the expected value.  Summing all these deviations seems like a good measure of the total deviation in all the bins but before doing the sum it probably would be a good idea to square each deviation so eliminate cancellations (plus and minu) of deviations from different bins.  So a measure of the goodness of fit could be

(∑_k)((((Subscript[O, k])-(Subscript[E, k])))^2)

where the sum is over all the bins.  However, this is not quite the standard measure of goodness of fit.

Since you expect the observed data to fit a Poisson distribution you also expect that the "width" of the probability distribution is σ = (Sqrt[(Subscript[E, k])]).  This comes from the idea that the standard deviation (or width) of the Poisson distribution is (Sqrt[N]).  (The standard deviation of most probability distributions is roughly (Sqrt[N]))  So we believe that IF THE MODEL DISTRIBUTION is a "good fit" of our experimental data then MOST (actually 67% of the data of average) of the observed data points (occurances) in the kth bin will be in the range

(Subscript[(E ), k]) ±  (Sqrt[(Subscript[E, k])])

Using equation (2) in each term of equation (1), we see that the expect size of each term of equation (1) is

((((Subscript[O, k])-(Subscript[E, k])))^2)≈ ((((Subscript[(E ), k]) ±  (Sqrt[(Subscript[E, k])])-(Subscript[E, k])))^2)≈ ((((Sqrt[(Subscript[E, k])])))^2)=(Subscript[E, k])

So if we use equation (1) is our measure of the "goodness of fit" then the expect size of the sum of equation (1) is roughly b (Subscript[E, k]) where b is the number of bin.  We would like to have a criteria for goodness of fit that is more general (so we can construct and use tables).  Thus we settle on the criteria

(χ^2) =(∑_k)(((((Subscript[O, k])-(Subscript[E, k])))^2)/((Subscript[E, k])))

It is a good idea to use (Subscript[E, k]) in the denominator of (Subscript[E, k]) since the theoretically obtained (Subscript[E, k]) can be small but never zero (The experimentally obtained (Subscript[O, k]) can be observed to be zero for some bins like bins k=7, 9, 10 in the data below).  Now each term of the sum in equation (4) has (approximate) size ONE provided the data is "good" since

(((((Subscript[O, k])-(Subscript[E, k])))^2)/((Subscript[E, k])))≈ (((((Sqrt[(Subscript[E, k])])))^2)/(Subscript[E, k]))≈1

and the sum in equation (4) namely (χ^2) is of order of b the total number of bins.  If the observational data (Subscript[O, k]) is NOT "good" (i.e. poorly or carelessly taken data) or the data is NOT well fit by the model theory, then we expect (χ^2) to be LARGER than b and if we have especially good data, (χ^2) will be less than b.

Workshop

Recall the background radiation where Ntotal=50 observations were made each 10sec long.  Suppose your experimentally observed data for o[k] the number of counts k (where k=0, 1, 2, 3, ...)  in 10 seconds is below:  (Use the sample data below today but when you do you lab write up, use the data you take experimentally.)

({{k, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, {O[k], 2, 15, 9, 8, 10, 3, 2, 0, 1, 0, 0}})

The number of observations is Num=50.

1.  Input the data in the table above using for example

o[0] = 2. ;
o[1] = 15. ;
o[2] = 9. ;

etc.

INPUTING DATA:  Make sure you select the cell bracket and hit the enter key so Mathematica knows the data.  If you just type in the data without hitting enter, Mathematica has not knowledge of your data.

After typing AND inputing the above data, it is a good idea to check to make sure Mathematica knows the data.  Type

?o

and Mathematica should report on what it knows (the specfic numerical values) regarding the variable o.

IMPORTANT:  Use a small letter o[k] instead of capitol O[k] since O is reserved by Mathematica for special use.  You could also use a combination of letters instead of just o.

Graph the observed number of occurences o[k] versus the number of counts. First use the Table function to define data via

ObsData=Table[{k,o[k]},{k, 0, 10}];

The data table is assigned to the variable ObsData for use in ListPlot below.  Check to see that ObsData is a list of numbers ready for graphing by typing and inputing

?ObsData

The actual graph is obtained using

ObsGraph=ListPlot[ObsData,AxesLabel -> {"k", "o[k]"},PlotStyle -> PointSize[0.02]]

Next use the data above for o[k] to calculate the experimentally observed probability Pobs[k]=o[k]/Num.  First input the formula for Po[k] below:

Pobs[k_]:=o[k]/Num

The "bar" indicated in red is important in defining functions using Mathematica.  Using a colon (also indicated in red) is also a good idea.

Graph Pobs[k] versus the number of counts k using the Mathematica Table function and ListPlot using what was done with the experimental occurences data as a template.  You can make the dots smaller to distinguish from them from the data graph.  Assign your graph the name.  Assign this graph the name PobsGraph for later use.

2.  Calculate the average number of counts λ observed experimentally using

λ = (∑_(k = 0)^10) k*Pobs[k]

The sum k is over the bin numbers k=0 through k=10.

3.  Let us assume your experimental data can be modeled reasonably well by the theoretical Poisson distribution.  Use λ the average number of counts observed experimentally and calculated in item #2 above to define the theoretical (or expected) Poisson distribution P[k]

P[k_]:=(((E^(-λ))*(λ^k))/(k!))

Graph P[k] versus the number of counts k.  First use the Table function to define TheoData and then use ListPlot to do the actual graphing:

TheoData=Table[{k, P[k]}, {k, 0, 10}];
TheoGraph=ListPlot[TheoData]

Assign the name TheoGraph to your graph of the Poisson distribution. Use the Mathematica functon Show to view the ExpGraph and TheoGraph together.  Distinguish one graph from the other.

4.  Calculate the expected number of occurences e[k] for have k counts based on the theoretical Poisson distribution P[k] using

e[k] = Num*P[k]

You cannot use the capitol E as in E[k] because Mathematica reserved E for the exponential function.

5.  Calculate the (χ^2) corresonding to your data set using

χsqr =(∑_(k = 0)^10)(((((o[k] - e[k]))^2))/(e[k]))

This is like equation (4) but the reserved Mathematica function O and E do not appear.

QUESTION:  
	(1) Is your (χ^2) less than the number of bins b where you have nonzero observations?  If (χ^2)≤ b then you have "good data" and what this really means is that your model and experimental data agree so you say you understand the experiment.  
	(2) If you (χ^2)> b  then theory and experiment do not agree.  There are two possible reasons for this disagreement.  Possibly your theory is not appropriate and you should look for another.  In this case, you say you do not yet understand your experiment.  Alternatively, perhaps you have taken "bad data" either by being careless in taking your data or in having some poorly controlled effect.  For example, if the Geiger scaler had some residual radioactive contamination from someone doing the half-life experiment, then the average number of backgrounds would be higher than expected.  (Alternatively the beta source could be too close to the Geiger counter producing counts in addition to background.)  If you suspect you have taken "bad data" you should go back to your experiment and either be more careful or try to control the experimental situation better.  If you  have been reasonably careful taking the experimental data and still (χ^2)  is large, then you should search for a new theory.

©  Rodney L. Varley (2007).