Elliot Tanis
April 19, 2006

(a) Let X  have a Poisson distribution with parameter . Show that the mean and variance are both equal to . Do this using sums and also using the moment-generating function. Also find the third and fourth moments.

 > t := 't':

 > f := PoissonPDF(lambda, x); # Poisson p.d.f.

 > mu := simplify(sum(x*f, x = 0 .. infinity));

 > var := simplify(sum((x - mu)^2*f, x = 0 .. infinity));

 > M := sum(exp(t*x)*f, x = 0 .. infinity); # moment-generating function

 > Mp := diff(M, t); mu := simplify(subs(t = 0, Mp));

 > Mpp := diff(M, t\$2); var := simplify(subs(t = 0, Mpp) - mu^2);

 > Mppp := diff(M, t\$3); ThirdM := simplify(subs(t = 0, Mppp));

 > Mpppp := diff(M, t\$4); FourthM := simplify(subs(t = 0, Mpppp));

(b) We know that the mean and variance of   X -bar are  and / n . Perhaps not as well known is that   the mean of the sample variance       is   , and the variance of      is    .   {Using   material from Cramer   (1946, p. 348),   it can be shown that the variance of      is      where      is the fourth moment about the mean.}

We first find the mean and variance of the sample mean using the moment-generating function of  X- bar .

 > n := 'n': M := t-> exp(lambda*n*(exp(t/n) - 1));

 > Mp := diff(M(t), t); mu := simplify(subs(t = 0, Mp));

 > Mpp := diff(M(t), t\$2); var := simplify(subs(t = 0, Mpp) - mu^2);

Now find the fourth moment of  X about     using a sum .

 > mu := simplify(sum(x*PoissonPDF(lambda, x), x = 0 .. infinity));

 > mu[2] := simplify(sum((x - lambda)^2*PoissonPDF(lambda, x), x = 0 .. infinity));

 > mu[4] := simplify(sum((x - lambda)^4*PoissonPDF(lambda, x), x = 0 .. infinity));

 > VarVar := (n^2/(n-1)^2)*((mu[4] - mu[2]^2)/n - 2*(mu[4] - 2*mu[2]^2)/n^2 + (mu[4] - 3*mu[2]^2)/n^3);

 > VarVar := simplify(%);

We now compare the above answers with empirical characteristics.

 > nn := 10; l := 5.9; # l = lambda, the mean

 > randomize(): for k from 1 to 400 do   X := PoissonS(l, nn):   xxbar[k] := Mean(X):   ssvar[k] := Variance(X): od: Xbar := [seq(xxbar[k], k = 1 .. 400)]: Svar := [seq(ssvar[k], k = 1 .. 400)]: xbar := evalf(Mean(Xbar)); mu := l; varxbar := evalf(Variance(Xbar));  # sample variance of x-bar's theovarxbar := evalf(l/nn);  # theoretical variance of x-bar svarbar := evalf(Mean(Svar));  # sample mean of svar's theovarmean := l;  # theoretical mean of sample variance svarsvar := evalf(Variance(Svar));  # sample variance of svar's theosvar := evalf(subs(lambda=l, n=nn, VarVar)); # theoretical svar variance

 > f := PoissonPDF(5.9,x);

Here is a check of the procedure  "PoissonS"  that simulates Poisson observations.

 > randomize(): X := PoissonS(5.9, 400): P1 := HistogramFill(X, -0.5 .. 16.5,17): P2 := ProbHistB(f, 0 .. 17): xtics := [seq(k, k = 0 .. 17)]: ytics := [seq(0.01*k, k = 1 .. 25)]: P3 := plot([[0,0],[0,0]], x = -0.5 .. 16.5, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): display({P1, P2, P3});

 >

For determining the half-lives of radioactive isotopes, it is important to know what the background radiation is in a given detector over a period of time. Data taken in a   - ray detection experiment over 300 one-second intervals yielded the following data:

 > X := [0, 2, 4, 6, 6, 1, 7, 4, 6, 1, 1, 2, 3, 6, 4, 2, 7, 4, 4, 2, 2, 5, 4, 4, 4, 1, 2, 4, 3, 2, 2, 5, 0, 3, 1, 1, 0, 0, 5, 2, 7, 1, 3, 3, 3, 2, 3, 1, 4, 1, 3, 5, 3, 5, 1, 3, 3, 0, 3, 2, 6, 1, 1, 4, 6, 3, 6, 4, 4, 2, 2, 4, 3, 3, 6, 1, 6, 2, 5, 0, 6, 3, 4, 3, 1, 1, 4, 6, 1, 5, 1, 1, 4, 1, 4, 1, 1, 1, 3, 3, 4, 3, 3, 2, 5, 2, 1, 3, 5, 3, 2, 7, 0, 4, 2, 3, 3, 5, 6, 1, 4, 2, 6, 4, 2, 0, 4, 4, 7, 3, 5, 2, 2, 3, 1, 3, 1, 3, 6, 5, 4, 8, 2, 2, 4, 2, 2, 1, 4, 7, 5, 2, 1, 1, 4, 1, 4, 3, 6, 2, 1, 1, 2, 2, 2, 2, 3, 5, 4, 3, 2, 2, 3, 3, 2, 4, 4, 3, 2, 2, 3, 6, 1, 1, 3, 3, 2, 1, 4, 5, 5, 1, 2, 3, 3, 1, 3, 7, 2, 5, 4, 2, 0, 6, 2, 3, 2, 3, 0, 4, 4, 5, 2, 5, 3, 0, 4, 6, 2, 2, 2, 2, 2, 5, 2, 2, 3, 4, 2, 3, 7, 1, 1, 7, 1, 3, 6, 0, 5, 3, 0, 0, 3, 3, 0, 2, 4, 3, 1, 2, 3, 3, 3, 4, 3, 2, 2, 7, 5, 3, 5, 1, 1, 2, 2, 6, 1, 3, 1, 4, 4, 2, 3, 4, 5, 1, 3, 4, 3, 1, 0, 3, 7, 4, 0, 5, 2, 5, 4, 4, 2, 2, 3, 2, 4, 6, 5, 5, 3, 4]:

 > xbar := evalf(Mean(X)); var := evalf(Variance(X));

 > f := PoissonPDF(xbar,x);

 > Min(X),Max(X);

 > P1 := HistogramFill(X, -0.5 .. 8.5,9): P2 := ProbHistB(f, 0 .. 9): xtics := [seq(k, k = 0 .. 9)]: ytics := [seq(0.02*k, k = 1 .. 14)]: P3 := plot([[0,0],[0,0]], x = -0.5 .. 8.5, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): display({P1, P2,P3});

The Normal Approximation of the Poisson Distribution

For what values of  does the   N( , ) distribution provide good approximations for the Poisson distribution with mean ?

 > xtics := [seq(2*k, k = 1 .. 20)]: ytics := [seq(0.01*k, k = 1 .. 15)]: for L from 7 to 25 do H := ProbHistFillY(PoissonPDF(L, x),  0 .. 40): N := plot(NormalPDF(L, L, x), x = 0 .. 40,color=blue, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): t := textplot([30, 0.12, `l =` ||L], font=[SYMBOL]): H||L := display({H, N, t}): od: display([seq(H||L,  L = 7 .. 25)], insequence = true);

Poisson Probability Histogram ,   N ( ,   ) p.d.f.

Click on the above figure and use the tool bar to animate it.

 >