I needed to compute the higher moments of a mixture distribution for a project I’m working on. I’m writing up the code here in case anyone else finds this useful. (And in case I’ll find it useful in the future.) I’ll include the central moments first. From there it’s easy to compute skewness and kurtosis.

Suppose *X* is a mixture of *n* random variables *X*_{i} with weights *w*_{i}, non-negative numbers adding to 1. Then the *j*th central moment of *X* is given by

where μ_{i} is the mean of *X*_{i}.

In my particular application, I’m interested in a mixture of normals and so the code below computes the moments for a mixture of normals. It could easily be modified for other distributions.

from scipy.misc import factorialk, comb def mixture_central_moment(mixture, moment): '''Compute the higher moments of a mixture of normal rvs. mixture is a list of (mu, sigma, weight) triples. moment is the central moment to compute.''' mix_mean = sum( w*m for (m, s, w) in mixture ) mixture_moment = 0.0 for triple in mixture: mu, sigma, weight = triple for k in range(moment+1): prod = comb(moment, k) * (mu-mix_mean)**(moment-k) prod *= weight*normal_central_moment(sigma, k) mixture_moment += prod return mixture_moment def normal_central_moment(sigma, moment): '''Central moments of a normal distribution''' if moment % 2 == 1: return 0.0 else: # If Z is a std normal and n is even, E(Z^n) == (n-1)!! # So E (sigma Z)^n = sigma^n (n-1)!! return sigma**moment * factorialk(moment-1, 2)

Once we have code for central moments, it’s simple to add code for computing skewness and kurtosis.

def mixture_skew(mixture): variance = mixture_central_moment(mixture, 2) third = mixture_central_moment(mixture, 3) return third / variance**(1.5) def mixture_kurtosis(mixture): variance = mixture_central_moment(mixture, 2) fourth = mixture_central_moment(mixture, 4) return fourth / variance**2 - 3.0

Here’s an example of how the code might be used.

# Test on a mixture of 30% Normal(-2, 1) and 70% Normal(1, 3) mixture = [(-2, 1, 0.3), (1, 3, 0.7)] print "Skewness = ", mixture_skew(mixture) print "Kurtosis = ", mixture_kurtosis(mixture)

**Related post**: General formula for normal moments

Thank you very much. Very useful for me.

Cool, do you think it is possible having e.g. mean, stddev skew and kurtosis, and to estimate the parameters of a double- or tripple normal mix from them?

No, because there are too many degrees of freedom. If you have a mixture of n normals, you have 3n – 1 degrees of freedom: n-1 weights (you lose one degree of freedom because the mixture weights must add to one) and a mean and standard deviation for each normal in the mix.

So you need at least 3n-1 pieces of data. Then the question becomes what those pieces should be. I doubt that it would be possible/practical to use moments.