- mean = ((src->mean * (double) (nr - 1))
- + dst->mean) / ((double) nr);
- S = ((src->S * (double) (nr - 1)) + dst->S) / ((double) nr);
+ double delta = src->mean - dst->mean;
+
+ mean = ((src->mean * src->samples) +
+ (dst->mean * dst->samples)) /
+ (dst->samples + src->samples);
+
+ S = src->S + dst->S + pow(delta, 2.0) *
+ (dst->samples * src->samples) /
+ (dst->samples + src->samples);