www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - Calculating mean and standard deviation with std.algorithm.reduce

reply Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
Hello all,

A little challenge that's been bothering me today.

The docs for std.algorithm give an illustration of its use to calculate mean
and 
standard deviation in a single pass:

       // Compute sum and sum of squares in one pass
       r = reduce!("a + b", "a + b * b")(tuple(0.0, 0.0), a);
       assert(r[0] == 35); // sum
       assert(r[1] == 233); // sum of squares
       // Compute average and standard deviation from the above
       auto avg = r[0] / a.length;
       auto stdev = sqrt(r[1] / a.length - avg * avg);

However, this formula for standard deviation is one that is well known for
being 
subject to potentially fatal rounding error.  Consider the following:

       float[] a = [10_000.0f, 10_001.0f, 10_002.0f];

       auto r = reduce!("a + b", "a + b * b")(tuple(0.0f, 0.0f), a);

       auto avg = r[0] / a.length;
       auto sd = sqrt(r[1] / a.length - avg * avg);
       writeln(avg, "\t", sd);

... which gives a value of 0 for standard deviation when compiled with LDC or 
GDC, and -nan with DMD, when the correct answer is 0.816497.

An improved online formula calculates together two quantities:

              {  x[1]    if k == 1
     M[k]  =  {
              {  M[k-1] + (x[k] - M[k-1]) / k    , if k > 1

              {  0      if k == 1
     Q[k]  =  {
              {  Q[k-1] + (k - 1) * ((x[k] - M[k-1]) ^^ 2) / k   ,  if k > 1

... which one can readily prove give you respectively the mean of x[1], ..., 
x[k] and the sum of squares of deviations from the mean; so that you can 
calculate sample variance as Q[n] / (n - 1), or standard variance as Q[n] / n 
with standard deviation being given by the square root of this.

In a reduce sense, you'd want to calculate the mean according to,

     reduce!"a + (b - a)/k"(0.0, x);

... where k represents the index count 1, 2, 3, ...  However, it's not evident 
to me how you could get reduce() to know this counting value.

Where calculating Q[k] is concerned, you need to know both the index value
_and_ 
the value of the corresponding M[k].  Again, it's not evident to me how you'd 
indicate to reduce() what is needed.

Can anyone offer advice on how to achieve these calculations using reduce?

Thanks & best wishes,

     -- Joe
Feb 13 2013
next sibling parent reply "jerro" <a a.com> writes:
 ... where k represents the index count 1, 2, 3, ...  However, 
 it's not evident to me how you could get reduce() to know this 
 counting value.


You would use zip and sequence to add indices to x, like this:

reduce!reducer(initial, zip(x, sequence!"n"))


 Where calculating Q[k] is concerned, you need to know both the 
 index value _and_ the value of the corresponding M[k].  Again, 
 it's not evident to me how you'd indicate to reduce() what is 
 needed.


I guess reduce would need to operate on tuples of M and Q, so you 
would have something like:

alias Tuple!(float, float) MQ;

MQ reducer(MQ mq, Tuple!(float, int) xi)
{
     return MQ(
         // compute new values or M and Q here
     );
}

And you would then call reduce like this:

reduce!reducer(MQ(x.front, 0), zip(x, sequence!"n"))

You could also use Tuple of M, Q and k instead of using zip and 
sequence. You would then pass MQk(x.front, 0, 0) as first 
argument to reduce (I'm assuming zero based indexing here) and 
simply compute the k component of the return value in reducer as 
mqk[2] + 1.
Feb 13 2013
parent reply "jerro" <a a.com> writes:
 reduce!reducer(MQ(x.front, 0), zip(x, sequence!"n"))


A small correction : you would need to use x.drop(1) instead of 
x, because the first element of x is only used to compute the 
initial value of 1. If you wanted k to have the same meaning as 
the one in your formula, you would need to use sequence!"n + 2" 
instead of sequence!"n".
Feb 13 2013
parent reply "bearophile" <bearophileHUGS lycos.com> writes:
jerro:


 reduce!reducer(MQ(x.front, 0), zip(x, sequence!"n"))


 A small correction : you would need to use x.drop(1) instead of 
 x, because the first element of x is only used to compute the 
 initial value of 1. If you wanted k to have the same meaning as 
 the one in your formula, you would need to use sequence!"n + 2" 
 instead of sequence!"n".


See enumerate():
http://d.puremagic.com/issues/show_bug.cgi?id=5550

I think with enumerate it becomes:
MQ(x.front, 0).enumerate(1).reduce!reducer()

Bye,
bearophile
Feb 13 2013
parent "jerro" <a a.com> writes:
 See enumerate():
 http://d.puremagic.com/issues/show_bug.cgi?id=5550


I like this enumerate() thing. Is there any particular reason why
it isn't in phobos, or is it just that no one has added it yet?


 I think with enumerate it becomes:
 MQ(x.front, 0).enumerate(1).reduce!reducer()


I think the argument to enumerate here should be 2 (and the x
range is missing, of course).

Another way to do this is:

MQk(x.front, 0, 2).reduce!reducer(x.drop(1))

Or using a lambda instead of reducer():

auto mqk = tuple(x.front, 0.to!double, 2).reduce!
      ((mqk, x) => tuple(
          mqk[0] + (x - mqk[0]) / mqk[2],
          mqk[1] + (mqk[2] - 1) * ((x - mqk[0]) ^^ 2) / mqk[2],
          mqk[2] + 1))
      (x.drop(1));
Feb 13 2013
prev sibling parent reply FG <home fgda.pl> writes:
On 2013-02-13 14:44, Joseph Rushton Wakeling wrote:

 The docs for std.algorithm give an illustration of its use to calculate mean
and
 standard deviation in a single pass: [...]
 However, this formula for standard deviation is one that is well known for
being
 subject to potentially fatal rounding error.


Typical thing with examples - they try to be terse and show off a mechanism
like 
reduce, without going into too much details and hence are unusable IRL.

You can use reduce and put the division and subtraction into the reduce itself 
to prevent overflows. You also won't end up with jaw-dropping tuples, sorry. :)

     float[] a = [10_000.0f, 10_001.0f, 10_002.0f];
     auto n = a.length;
     auto avg = reduce!((a, b) => a + b / n)(0.0f, a);
     auto var = reduce!((a, b) => a + pow(b - avg, 2) / n)(0.0f, a);
     auto sd = sqrt(var);
     writeln(avg, "\t", sd);

Output: 10001   0.816497
Feb 13 2013
next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 02/13/2013 03:48 PM, FG wrote:

 You can use reduce and put the division and subtraction into the reduce itself
 to prevent overflows. You also won't end up with jaw-dropping tuples, sorry. :)

      float[] a = [10_000.0f, 10_001.0f, 10_002.0f];
      auto n = a.length;
      auto avg = reduce!((a, b) => a + b / n)(0.0f, a);
      auto var = reduce!((a, b) => a + pow(b - avg, 2) / n)(0.0f, a);
      auto sd = sqrt(var);
      writeln(avg, "\t", sd);

 Output: 10001   0.816497


Nice :-)

The thing that's desirable about the example given in the docs is that it's a 
one-pass approach to getting standard deviation (or rather, mean and standard 
deviation together).  Two reduce commands as you suggest is fine for small
data, 
but if you have something really large you'd rather compute it in one pass. 
You 
might also have input that is not [easily] repeatable, so you _have_ to do it 
online -- imagine that your input range is e.g. a file byLine (you'd prefer to 
avoid parsing that twice), or input from a tape (OK, OK, old fashioned) or 
perhaps data from an external source.
Feb 13 2013
prev sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 02/13/2013 03:48 PM, FG wrote:

 Typical thing with examples - they try to be terse and show off a mechanism
like
 reduce, without going into too much details and hence are unusable IRL.


My favourite -- in the tutorial for a really serious piece of scientific code 
written in C:

      int n = rand() % 10;
Feb 13 2013