// from https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm

srecord Welford(long count, dbl mean, dbl m2) {

  // For a new value newValue, compute the new count, new mean, the new M2.

  // mean accumulates the mean of the entire dataset

  // M2 aggregates the squared distance from the mean

  // count aggregates the number of samples seen so far

  void add(double newValue) {

    count++;

    dbl delta = newValue - mean;

    mean += delta / count;

    dbl delta2 = newValue - mean;

    m2 += delta * delta2;

  }

  void remove(double valueToRemove) {

    guarantee(count > 0);

    /* 

      Let's reverse the formulas from add():

        delta = newValue - mean1

        mean2 = mean1 + delta / count2

        => mean2 = mean1+(newValue-mean1)/count2

                 = mean1 + newValue/count2 - mean1/count2

        => mean2 = mean1*(1-1/count2) + newValue/count2

      Solve for mean1. (| -newValue/count2)

        => mean2-newValue/count2 = mean1*(1-1/count2)  | :(1-1/count2)

        => (mean2-newValue/count2)/(1-1/count2) = mean1.

      "Count" is easy to reverse: count1 = count2-1, obviously.

      Now, about restoring m2.

        delta2 = newValue - mean2

        m2_2 = m2_1 + delta * delta2

        => m2_2 - delta * delta2 = m2_1

      Well, that was easy!

    */

    double newValue = valueToRemove;

    double delta2 = newValue-mean;

    mean = (mean-newValue/count)/(1-1.0/count);

    double delta = newValue-mean;

    m2 -= delta*delta2;

    --count;

  }

  bool isEmpty() { ret count == 0; }

  dbl average aka avg aka mean() { ret count == 0 ? Double.NaN : mean; }

  dbl variance() { ret count == 0 ? Double.NaN : m2/count; }

  dbl sampleVariance() { ret count < 2 ? Double.NaN : m2/(count-1); }

  dbl standardDeviation() { ret sqrt(variance()); }

  toString {

    S s = nValues(count);

    if (!isEmpty())

      s += ", avg: " + avg() + ", std dev: " + standardDeviation();

    ret s;

  }

  AverageAndStandardDeviation get() {

    ret new AverageAndStandardDeviation(avg(), standardDeviation());

  }

}