diff --git a/README.md b/README.md index 7f94108e463ffd135eee53a27e4b296cb1d98a85..52e8fa69ddfb8fe5efc6a90ab6e90bdeaa60db94 100644 --- a/README.md +++ b/README.md @@ -566,20 +566,20 @@ the uncertainties are the square roots of the expected number of events. If this case, we have to compute the $`\chi^2`$ manually. The following example shows both, the usage of `chisquare` and the manual computation. The print statements for both methods produce the same output. Please note that -our model has five degrees of freedom since we have five free parameters. The +our model has five free parameters. The total number of degrees of freedom is the number of data points reduced by the -number of degrees of freedom of our model. +of parameters of our model. <!-- append decay.py --> ```python -# Degrees of freedom of our model. -dof = 5 +# Delta degrees of freedom +ddof = 5 print("chi^2 from scipy:") # Call chisquare and print. -chi2, p = scipy.stats.chisquare(count, fit_count, ddof=dof) -print(" chi2 / dof = %g / %d" % (chi2, len(count) - dof)) +chi2, p = scipy.stats.chisquare(count, fit_count, ddof=ddof - 1) +print(" chi2 / dof = %g / %d" % (chi2, len(count) - ddof)) print(" p-value = %g" % p) print() # print blank line @@ -591,17 +591,17 @@ uncertainty = np.sqrt(fit_count) chi2 = ((count - fit_count)**2 / uncertainty**2).sum() # Calculate p-value and print -p = scipy.stats.distributions.chi2.sf(chi2, len(count) - 1 - dof) -print(" chi2 / dof = %g / %d" % (chi2, len(count) - dof)) +p = scipy.stats.distributions.chi2.sf(chi2, len(count) - ddof) +print(" chi2 / dof = %g / %d" % (chi2, len(count) - ddof)) print(" p-value = %g" % p) - ``` -The additional `-1` stems from the fact that we have a Poisson distribution -for which the uncertainties are derived from the number of events. The -test is not sensitive to the normalization, and therefore, we loose one degree -of freedom, cf. [Data Analysis in High Energy -Physics](https://www.wiley.com/en-us/Data+Analysis+in+High+Energy+Physics%3A+A+Practical+Guide+to+Statistical+Methods-p-9783527410583) by Behnke et at. +The additional `- 1` in the call of `chisquare` is necessary because the method +implicitly assumes that the total number of events is fixed. The method +therefore reduces the number of degrees of freedom by one. However, in our case, +the number of decays is not fixed and we have an additional degree +of freedom compared to what `chisquare` assumes. The `- 1` in the `ddof` +corrects for this. If you run the `decay.py` you should see the following fit results. <!-- console_output --> @@ -616,11 +616,11 @@ Optimal parameters: chi^2 from scipy: chi2 / dof = 120.577 / 123 - p-value = 0.519418 + p-value = 0.544943 Manual chi^2 test: chi2 / dof = 120.577 / 123 - p-value = 0.519418 + p-value = 0.544943 ``` Congratulations! You have mastered the first steps to analyze experimental data