Higher moments of normal distribution

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 Xn by taking the nth 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

Tagged with: ,
Posted in Python
5 comments on “Higher moments of normal distribution
  1. Simon H says:

    Nice snippet!

    FWIW Wolfram Alpha will also provide them now.

  2. GlennF says:

    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.

  3. GlennF says:

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

  4. Jan Van lent says:

    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} $

  5. Jan Van lent says:

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