# Formulae for combining statistical moments

I am writing code to calculate statistical moments (mean, variance, skewness, kurtosis) for large samples of data and have the requirement of needing to be able to calculate moments for subsections of the sample (in parallel), then combine/merge them to give the moment for the sample as a whole.

For example:

$S = \lbrace 1.0, 1.2, 2.0, 1.7, 3.4, 0.9 \rbrace$

$A = \lbrace 1.0, 1.2, 2.0 \rbrace$ and $B = \lbrace 1.7, 3.4, 0.9 \rbrace$

So $A \cup B = S$

I need to calculate the statistics/moments for $A$ and $B$, then combine them to give the statistics/moments for $S$

Count is simple: $n_S = n_A + n_B$

Mean is not much worse: $\mu_S = (n_A\mu_A + n_B\mu_B) / n_S$

Variance is a little less pretty: $\sigma_S = [n_A\sigma_A + n_B\sigma_B + (\frac{n_An_A}{n_A+n_B})(\mu_A - \mu_B)^2] / n_S$

But now I'm struggling for skewness and, in particular, kurtosis. I have all 'lesser' moments for each of the subsections of the data available and have some idea of the direction I'm heading, but am really struggling to derive the formulae needed.

Has anybody derived these formulae before? Could anyone point me in the right direction? These may be simple/obvious things to any with anyone with a decent amount of statistical knowledge, unfortunately that's something I completely lack...

-

I happened to solve exactly this problem at my previous job.

Given samples of size $n_A$ and $n_B$ with means $\mu_A$ and $\mu_B$, and you want to calculate the mean, variance etc for the combined set $X=A\cup B$. A pivotal quantity is the difference in means

$$\delta = \mu_B - \mu_A$$

This already appears in your formula for variance. You could re-write your formula for the mean to include it as well, although I won't. I will, however, re-write your formulas to work with extensive terms (sums, sums of squares) rather than intensive terms (means, variances)

$$S_X = n_X\mu_X = n_A\mu_A + n_B\mu_B = S_A + S_B$$

$$S^2_X = n_X \sigma_X^2 = n_A\sigma_A^2 + n_B\sigma_B^2 + \frac{n_A n_B}{n_X} \delta^2 = S^2_A + S^2_B + \frac{n_A n_B}{n_X} \delta^2$$

Note that $S^j_X$ is the sum of differences from the mean, to the power $j$.

The formula for the sum of third powers, $S^3_X$, is

$$S^3_X = S^3_A + S^3_B + \frac{n_A n_B (n_A-n_B)}{n^2_X} \delta^3 + 3 \frac{n_A S^2_B - n_B S^2_A}{n_X} \delta$$

and for the sum of fourth powers

$$S^4_X = S^4_A + S^4_B + \frac{n_A n_B (n_A^2 - n_A n_B + n_B^2)}{n^3_X} \delta^4 + 6\frac{n^2_A S^2_B + n^2_B S^2_A}{n^2_X} \delta^2 + 4\frac{n_A S^3_B - n_B S^3_A}{n_X} \delta$$

Once you have these quantities, you can calculate the quantities you're interested in:

$$\mu_X = \frac{S_X}{n_X}$$

$$\sigma^2_X = \frac{S^2_X}{n_X}$$

$$s_X = \frac{\sqrt{n_X}S^3_X}{(S^2_X)^{3/2}}$$

$$\kappa_X = \frac{n_X S^4_X}{(S^2_X)^2}$$

Needless to say, you should write unit tests that compare the output from these formulas to the ones computed in the 'traditional' way to make sure that you (or I) haven't made a mistake somewhere :)

-
Hey Chris, this is a great answer, but not necessarily what I need in my case (which I should have clarified in my question) due to massive cancellation. The formulae need to be in terms of intensive terms and avoid enormous quantities (such as $S^4_X$) to avoid loss of precision. Essentially I need $S^4_X$ expressed in terms of central moments / intensive quantities... –  Simon Cowen Oct 19 '12 at 11:48
Hi @Simon. That's no problem - each of the equations can be divided by $n_X$ to get the corresponding intensive statistic, and then you can use the four equations at the end to remove all references to $S_X$, $S^2_X$ etc in favour of $\mu_X$, $\sigma_X^2$ etc. However, this isn't where you should be worrying about cancellation errors. The quantity $\delta$ is constructed by subtracting two similar-sized floating point numbers, so you expect some loss of precision there (how much depends on how large the means are, and how much they differ between samples). You could re-write all the... –  Chris Taylor Oct 19 '12 at 11:58
...formulas to expand out $\delta$ as a sum of terms involving $\mu_A$ instead, and cancel as many terms manually as you can (in fact you should use a tool like Mathematica or Wolfram Alpha to do this instead of doing it manually). You can then get a result that is less prone to loss of precision, at the expense of being more complicated to code up, and probably a bit slower to calculate. I gave you the mathematically prettiest version, rather than the version most suited to implementation, but it's not too much work to make the necessary changes. –  Chris Taylor Oct 19 '12 at 12:01
Saying all that, I have avoided the most pernicious problems with cancellation error when calculating statistical moments, by always working with terms that are differences from the mean. The biggest source of numerical error in calculating variance comes from using formula like $(1/n)\sum x_i^2 - \bar{x}^2$ instead of $(1/n)\sum (x_i-\bar{x})^2$. You might find something of value in this paper by Tony Finch, which goes through some of these considerations in detail. –  Chris Taylor Oct 19 '12 at 12:08
Chris, thanks for all this, you've certainly saved me some hours of deriving this myself. Thanks for the link to the paper too, I'm already using incremental algorithms but some of the other items may well come in handy later on. –  Simon Cowen Oct 22 '12 at 8:20