Running stats

running_stats.py

class RunningStats:
    """Computes running mean and standard deviation
    Adapted from:
        *
        <http://stackoverflow.com/questions/1174984/how-to-efficiently-\
calculate-a-running-standard-deviation>
        * <http://mathcentral.uregina.ca/QQ/database/QQ.09.02/carlos1.html>
    """

    def __init__(self):
        self.n = 0.

    def clear(self):
        self.n = 0.

    def push(self, x, per_dim=True):
        x = numpy.array(x).copy().astype('float16')
        # process input
        if per_dim:
            self.update_params(x)
        else:
            for el in x.flatten():
                self.update_params(el)

    def update_params(self, x):
        self.n += 1
        if self.n == 1:
            self.m = x
            self.s = 0.
        else:
            prev_m = self.m.copy()
            self.m += (x - self.m) / self.n
            self.s += (x - prev_m) * (x - self.m)

    def mean(self):
        return self.m if self.n else 0.0

    def variance(self):
        return self.s / (self.n) if self.n else 0.0

    def std(self):
        return numpy.sqrt(self.variance())
        
def test_running_stats():
    from numpy.testing import assert_almost_equal as almost_equal
    from numpy.random import randint
    arr = numpy.random.randn((30*8*12)).reshape((30, 8, 12))
    varsize_arr = []
    for el in arr:
        s = el.shape
        varsize_arr.append(el[:randint(s[0]-2)+2, :randint(s[1]-2)+2])

    # test per dimension statistics
    perdim_runner = RunningStats()
    for i, el in enumerate(arr, 1):
        perdim_runner.push(el)
        if i == 1:
            # arr[:i] has no axis 0
            continue
        almost_equal(arr[:i].mean(axis=0), perdim_runner.mean())
        almost_equal(arr[:i].std(axis=0), perdim_runner.std())

    # test single number statistics
    runner = RunningStats()
    for i, el in enumerate(varsize_arr, 1):
        runner.push(el, False)
        cum_arr = []
        for im in varsize_arr[:i]:
            cum_arr = numpy.concatenate([im.flatten(), cum_arr])
        almost_equal(numpy.array(cum_arr).mean(), runner.mean())
        almost_equal(numpy.array(cum_arr).std(), runner.std())