$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 2008-2011 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is free software; you can redistribute it and/or modify    */
/* it under the terms of the GNU General Public License as published by   */
/* the Free Software Foundation, either version 3 of the License, or      */
/* (at your option) any later version.                                    */

/* QUB Express is distributed in the hope that it will be useful,         */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of         */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          */
/* GNU General Public License for more details.                           */

/* You should have received a copy of the GNU General Public License,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

#include <math.h>
#include "running_mean.h"

/*
Here is a technique to estimate mean and std.dev of an unbounded sequence::
 The basic equation you need for a single pass update is:
 new S2 = old S2 + (x - old_M).(x - new_M)
 where S2 = sum of squares of deviations from the mean M.
 
 It is implemented as follows:
 Start: n = 0, M = 0, S2 = 0
 Then for each observation x:
 dev = x - M
 n = n + 1
 M = M + dev/n
 S2 = S2 + dev.(x - M)
 
 When you want to calculate the sample std. devn. :
 std.dev. = sqrt( S2 / (n-1) )
 --
 Alan Miller, Retired Scientist (Statistician)
 CSIRO Mathematical & Information Sciences
 a...@vic.cmis.csiro.au, mille...@ozemail.com.au
 http://www.ozemail.com.au/~milleraj
*/


RunningMean::RunningMean()
  : n(0), mean(0.0), s2(0.0)
{}

void RunningMean::add(double x)
{
  double dev = x - mean;
  ++n;
  mean += dev/n;
  s2 += dev*(x - mean);
}

double RunningMean::std()
{
  if ( n == 0 )
    return 0.0;
  if ( n == 1 )
    return mean;
  return sqrt(s2 / (n-1));
}