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...
 A: 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 :)
A: Ten years late here, but I wanted to offer a more general solution I recently came to for this problem. It turns out that the combined raw moments are simply weighted averages of the raw moments of each component set (A, B, C, etc. if you wanted to do this for more than two sets).
$$
\langle X^k \rangle = \frac{1}{n_{\rm tot}}\sum\limits_{i \in {A,B,...}}n_i \langle X^k_i \rangle
$$
This is why it is trivial to calculate the combined mean, since it is the first raw moment $\mu \equiv \langle X \rangle$ (note the brackets just symbolize the expectation value). For higher central/standardized moments like variance, skewness, and kurtosis, you just need to apply a binomial expansion to the raw moments:
$$
{\rm Var}(X) \equiv \sigma^2 = \langle (X - \mu)^2 \rangle \\ = \langle X^2 \rangle - \mu^2
$$
$$
{\rm Skew}(X) = \frac{1}{\sigma^3} \langle (X - \mu)^3 \rangle \\ = \frac{1}{\sigma^3} (\langle X^3 \rangle - 3\mu \langle X^2 \rangle + 2 \mu^3)
$$
$$
{\rm Kurtosis}(X) = \frac{1}{\sigma^4} \langle (X - \mu)^4 \rangle \\ = \frac{1}{\sigma^4} (\langle X^4 \rangle - 4\mu \langle X^3 \rangle + 6\mu^2 \langle X^2 \rangle - 3\mu^4)
$$
To be fully general, the combined n$^{\rm th}$ central moment is:
$$
\langle (X - \mu)^n \rangle = \sum\limits_{k=0}^n {n\choose k} (-1)^{n-k} \langle X^k \rangle \mu^{n-k}
$$
