>    read `C:\\Tanis-Hogg\\Maple Examples\\stat.m`:

Elliot Tanis
April 19, 2006

Poisson Examples

(a) Let X  have a Poisson distribution with parameter lambda . Show that the mean and variance are both equal to lambda . 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.

f := lambda^x*exp(-lambda)/x!

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

mu := lambda

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

var := lambda

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

M := exp(-lambda+exp(t)*lambda)

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

Mp := exp(t)*lambda*exp(-lambda+exp(t)*lambda)

mu := lambda

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

Mpp := exp(t)*lambda*exp(-lambda+exp(t)*lambda)+exp(t)^2*lambda^2*exp(-lambda+exp(t)*lambda)

var := lambda

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

Mppp := exp(t)*lambda*exp(-lambda+exp(t)*lambda)+3*exp(t)^2*lambda^2*exp(-lambda+exp(t)*lambda)+exp(t)^3*lambda^3*exp(-lambda+exp(t)*lambda)

ThirdM := lambda+3*lambda^2+lambda^3

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

Mpppp := exp(t)*lambda*exp(-lambda+exp(t)*lambda)+7*exp(t)^2*lambda^2*exp(-lambda+exp(t)*lambda)+6*exp(t)^3*lambda^3*exp(-lambda+exp(t)*lambda)+exp(t)^4*lambda^4*exp(-lambda+exp(t)*lambda)

FourthM := lambda+7*lambda^2+6*lambda^3+lambda^4

(b) We know that the mean and variance of   X -bar are lambda  and lambda / n . Perhaps not as well known is that   the mean of the sample variance    S^2    is   lambda , and the variance of    S^2   is   lambda*(2*lambda*n+n-1)/(n-1)/n  .   {Using   material from Cramer   (1946, p. 348),   it can be shown that the variance of    S^2   is   (mu[4]-(n-3)/(n-1)*sigma^4)/n    where    mu[4]   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));

M := proc (t) options operator, arrow; exp(lambda*n*(exp(t/n)-1)) end proc

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

Mp := lambda*exp(t/n)*exp(lambda*n*(exp(t/n)-1))

mu := lambda

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

Mpp := lambda/n*exp(t/n)*exp(lambda*n*(exp(t/n)-1))+lambda^2*exp(t/n)^2*exp(lambda*n*(exp(t/n)-1))

var := lambda/n

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

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

mu := lambda

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

mu[2] := lambda

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

mu[4] := 3*lambda^2+lambda

>    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 := n^2/(n-1)^2*((2*lambda^2+lambda)/n-2*(lambda^2+lambda)/n^2+lambda/n^3)

>    VarVar := simplify(%);

VarVar := (2*n*lambda+n-1)/n*lambda/(n-1)

We now compare the above answers with empirical characteristics.

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

nn := 10

l := 5.9

>    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

xbar := 5.846250000

mu := 5.9

varxbar := .5729683584

theovarxbar := .5900000000

svarbar := 5.715638889

theovarmean := 5.9

svarsvar := 8.237130040

theosvar := 8.325555555

>    f := PoissonPDF(5.9,x);

f := .2739444819e-2*5.9^x/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});

[Maple Plot]

>   

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   gamma - 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));

xbar := 3.030000000

var := 3.193076923

>    f := PoissonPDF(xbar,x);

f := .4831563813e-1*3.030000000^x/x!

>    Min(X),Max(X);

0, 8

>    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});

[Maple Plot]

The Normal Approximation of the Poisson Distribution

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

>    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);

[Maple Plot]

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

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

>   

Return to Menu