I want to use Welford's method to compute a running variance and mean. I came across this implementation of Welford's method in Python. However, when testing to double-check that it results in the same output as the standard Numpy implementation of calculating variance, I do find that there is a difference in output.
Running the following code (using the python module unittest) shows that the two give different results (even after testing many times):
random_sample = np.random.normal(0, 1, 100)
std = np.var(random_sample, dtype=np.longdouble)
mean = np.mean(random_sample, dtype=np.longdouble)
welford = Welford()
welford.add_all(random_sample)
self.assertAlmostEqual(mean, welford.mean)
self.assertAlmostEqual(var, welford.var_s)
>> AssertionError: 1.1782075496578717837 != 1.1901086360180526 within 7 places (0.011901086360180828804 difference)
Interestingly, there is only a difference in the variance, not the mean.
For my purposes, a 0.012 difference is significant enough that it could affect my results.
Why would there be such a difference? Could this be due to compounding floating point errors? If so, would my best bet be to rewrite the package to use the Decimal class?
By default, np.var computes the so-called "population variance", in which the number of degrees of freedom is equal to the number of elements in the array.
wellford.var_s is the sample variance, in which the number of degrees of freedom is the number of elements in the array minus one.
To eliminate the discrepancy, pass ddof=1 to np.var:
import numpy as np
from welford import Welford
random_sample = np.random.normal(0, 1, 100)
var = np.var(random_sample, dtype=np.longdouble, ddof=1)
welford = Welford()
welford.add_all(random_sample)
np.testing.assert_allclose(var, welford.var_s)
Alternatively, if it is appropriate to use the population variance in your application, use welford.var_p.
var = np.var(random_sample, dtype=np.longdouble)
np.testing.assert_allclose(var, welford.var_p)
For a description of the difference between the two, see the development version of the np.var documentation.