I have the following method of calculating the mean and standard deviation of a stream of numbers
import numpy as np
class RunningStatisticsVar:
def __init__(self, ddof=0):
self.mean = 0
self.var = 0
self.std = 0
self._n = 0
self._s = 0
self._ddof = ddof
def update(self, values):
values = np.array(values, ndmin=1)
n = len(values)
self._n += n
old_mean = self.mean
delta = values - self.mean
self.mean += (delta / self._n).sum()
self._s += (delta * (values - self.mean)).sum()
self.var = self._s / (self._n - self._ddof) if self._n > self._ddof else 0
self.std = np.sqrt(self.var)
def update_single(self, value):
self._n += 1
old_mean = self.mean
self.mean += (value - old_mean) / self._n
self._s += (value - old_mean) * (value - self.mean)
self.var = self._s / (self._n - self._ddof) if self._n > self._ddof else 0
self.std = np.sqrt(self.var)
def __str__(self):
if self.std:
return f"(\u03BC \u00B1 \u03C3): {self.mean} \u00B1 {self.std}"
else:
return f"{self.mean}"
Since I'd like to use this class in general applications, I want to make it as robust as possible. I've done my best to make the code numerically stable by trying to always place the summations and multiplications (i.e. the operations which can blow up numbers) as late as possible.
I have the below code, which I wrote to test the correctness of my implementation of the batch update in update from the update_single, which I knew to be correct:
samples = np.random.uniform(0, 100, [100000]) #
s1 = RunningStatisticsVar()
s2 = RunningStatisticsVar()
s1.update(samples)
print(s1)
for i in samples:
s2.update_single(i)
print(s2)
To try and test its numerical stability, I tried using some extremely small and some extremely large numbers by replacing the #-marked line with:
samples = np.hstack((np.random.uniform(0, 1e-20, [100000]), np.random.uniform(1e20, 1e21, [100000]), np.random.uniform(0, 1e-20, [100000])))
This code prints out:
(μ ± σ): 1.8363120959956487e+20 ± 3.000290391225753e+20
(μ ± σ): 1.836312095995585e+20 ± 3.000290391225752e+20
As you can see, the numbers are very much similar, but not quite the same. Is my code numerically stable and is there a better way to test it than what I have above?
