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

• Joseph Rushton Wakeling (41/41) Feb 13 2013 Hello all,
• jerro (18/25) Feb 13 2013 I guess reduce would need to operate on tuples of M and Q, so you
• jerro (5/6) Feb 13 2013 A small correction : you would need to use x.drop(1) instead of
• bearophile (7/13) Feb 13 2013 See enumerate():
• jerro (13/17) Feb 13 2013 I think the argument to enumerate here should be 2 (and the x
• FG (12/16) Feb 13 2013 Typical thing with examples - they try to be terse and show off a mechan...
```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 == 35); // sum
assert(r == 233); // sum of squares
// Compute average and standard deviation from the above
auto avg = r / a.length;
auto stdev = sqrt(r / 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 / a.length;
auto sd = sqrt(r / 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    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, ...,
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
``` ... 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 + 1.
```
Feb 13 2013
``` 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"
```
Feb 13 2013
```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"

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
``` 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 + (x - mqk) / mqk,
mqk + (mqk - 1) * ((x - mqk) ^^ 2) / mqk,
mqk + 1))
(x.drop(1));
```
Feb 13 2013
```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
```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    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