Sometimes a little bit of Python beats a Google search.

Last week I needed to look up the moments of a normal distribution. The first two moments are common knowledge, the next two are easy to find, but I wasn’t able to find the higher moments.

Here is a little Sage code that produces a table of moments for the normal distribution. (Sage is a Python-based mathematical computing environment.) The code computes the expected value of *X*^{n} by taking the *n*th derivative of the moment generating function and setting its argument to zero.

var('m, s, t') mgf(t) = exp(m*t + t^2*s^2/2) for i in range(1, 11): derivative(mgf, t, i).subs(t=0)

Here’s the output.

m m^2 + s^2 m^3 + 3*m*s^2 m^4 + 6*m^2*s^2 + 3*s^4 m^5 + 10*m^3*s^2 + 15*m*s^4 m^6 + 15*m^4*s^2 + 45*m^2*s^4 + 15*s^6 m^7 + 21*m^5*s^2 + 105*m^3*s^4 + 105*m*s^6 m^8 + 28*m^6*s^2 + 210*m^4*s^4 + 420*m^2*s^6 + 105*s^8 m^9 + 36*m^7*s^2 + 378*m^5*s^4 + 1260*m^3*s^6 + 945*m*s^8 m^10 + 45*m^8*s^2 + 630*m^6*s^4 + 3150*m^4*s^6 + 4725*m^2*s^8 + 945*s^10

**Update**: Here is a formula for the higher moments of the normal.

**Related post**: Sage Beginner’s Guide

Nice snippet!

FWIW Wolfram Alpha will also provide them now.

It’s kind of fun (and not too bad) to work it out in general. For the central case,

the expectation of X^n is zero if n is odd, and if n it even it is

(s^p p!)/((2^(p/2))*((p/2)!)).

For the non-central case, basically just use the binomial theorem. Sum p from 0 to n (over even values only) of:

Binom(n,p)*(s^p p!)/((2^(p/2))*((p/2)!))* m^(n-p)

Apologies if my hideous formating may make it look more complicated than it is.

Oops… substitute ‘n’ for ‘p’ in the central case expresssion.

If you don’t know the moment generating function you can also use Sage to evaluate the integrals for the moments directly. The

`integral`

function uses Maxima behind the scenes.var('m s t x')

assume(s > 0)

f = exp(-1/2*((x-m)/s)^2)

mom = lambda i: integral(x^i*f, (x, -oo, oo))

for i in range(1, 10+1):

print expand(mom(i)/mom(0))

With the same results as above. For LaTeX output:

print ",n".join("$latex {0} $".format(latex(expand(mom(i)/mom(0)))) for i in range(1, 10+1))

$latex 1 $,

$latex m $,

$latex m^{2} + s^{2} $,

$latex m^{3} + 3 , m s^{2} $,

$latex m^{4} + 6 , m^{2} s^{2} + 3 , s^{4} $,

$latex m^{5} + 10 , m^{3} s^{2} + 15 , m s^{4} $,

$latex m^{6} + 15 , m^{4} s^{2} + 45 , m^{2} s^{4} + 15 , s^{6} $,

$latex m^{7} + 21 , m^{5} s^{2} + 105 , m^{3} s^{4} + 105 , m s^{6}

$,

$latex m^{8} + 28 , m^{6} s^{2} + 210 , m^{4} s^{4} + 420 , m^{2}

s^{6} + 105 , s^{8} $,

$latex m^{9} + 36 , m^{7} s^{2} + 378 , m^{5} s^{4} + 1260 , m^{3}

s^{6} + 945 , m s^{8} $,

$latex m^{10} + 45 , m^{8} s^{2} + 630 , m^{6} s^{4} + 3150 , m^{4}

s^{6} + 4725 , m^{2} s^{8} + 945 , s^{10} $

Oops. I though WordPress would typeset the LaTeX output. Also the print statement in the loop should be indented.