From 29aa458b998ac4ddc71917e5dc8702e5dbee0096 Mon Sep 17 00:00:00 2001
From: Frank Sauerburger <frank@sauerburger.com>
Date: Thu, 19 Jul 2018 17:05:41 +0200
Subject: [PATCH] Correct -1 in chi2 term
---
README.md | 32 ++++++++++++++++----------------
1 file changed, 16 insertions(+), 16 deletions(-)
diff --git a/README.md b/README.md
index 7f94108..52e8fa6 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
--
GitLab