Typo in ci.confidence_interval
Closed this issue · 2 comments
Hello, I noticed a discrepancy in the the confidence intervals determined using allantools vs. those determined with matlab and checked against the NIST reference noted. It looks like there is a typo in ci.confidence_interval where the chi2_h and the chi2_l got swapped between var_l and var_h. The code snippet is below where the issue is in the last three lines of code
def confidence_interval(dev, edf, ci=one_sigma_ci):
""" returns confidence interval (dev_min, dev_max)
for a given deviation dev, equivalent degrees of freedom edf,
and degree of confidence ci.
Parameters
----------
dev: float
Mean value (e.g. adev) around which we produce the confidence interval
edf: float
Equivalent degrees of freedon
ci: float, defaults to scipy.special.erf(1/math.sqrt(2))
for 1-sigma standard error set
ci = scipy.special.erf(1/math.sqrt(2))
= 0.68268949213708585
Returns
-------
(dev_min, dev_max): (float, float)
Confidence interval
"""
ci_l = min(np.abs(ci), np.abs((ci-1))) / 2
ci_h = 1 - ci_l
# function from scipy, works OK, but scipy is large and slow to build
chi2_l = scipy.stats.chi2.ppf(ci_l, edf)
chi2_h = scipy.stats.chi2.ppf(ci_h, edf)
variance = dev*dev
var_l = float(edf) * variance / chi2_h # NIST SP1065 eqn (45)
var_h = float(edf) * variance / chi2_l
return (np.sqrt(var_l), np.sqrt(var_h))
In the output, this error is causing the confidence intervals to be swapped where the upper confidence interval is smaller than the lower confidence interval, where I would expect that due to the chi squared distribution being asymmetric when the mean is near zero, the lower confidence interval should be the smaller one.
the definition of ci_l with min() might be the problem here?
however the example code with one_sigma_ci seems to work: https://github.com/aewallin/allantools/blob/master/examples/ci_demo.py
could you please provide a code example that shows this problem?
My apologies, in the process of making the example code I realized it was an error in my understanding and the original was correct. Thank you for your discussion, I will close the comment!