1 . The Poisson distribution can be used to approximate binomial probabilities. In particular, if we take the limit of a b(n , p ) moment-generating function as proc (n) options operator, arrow; infinity end proc... such that np = lambda , a constant, then the limiting moment-generating function is that for the Poisson distribution with mean lambda .

> t := 't':
lambda := 'lambda':
Mp := exp(lambda*(exp(t)-1)); # Poisson mgf, lambda=2

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

> p := 'p':
n := 'n':
Mb := (1-p + p*exp(t))^n; # b(n,p) mgf

Mb := (1-p+p*exp(t))^n

Let lambda = 5 and graph the binomial moment-generating functions for the pairs: n = 10, p = 1/2; n = 20, p = 1/4, n = 50, p = 1/10.

> xtics := [seq(-0.40 + k*0.05, k = 0 .. 16)]:
ytics := [k$k=1..11]:

> lambda := 5:
Pp := plot(Mp, t = -.4 .. .4, xtickmarks=xtics, ytickmarks=ytics, labels =[``,``], color=black, thickness=2):
txtp := textplot([0.1,11,`Poisson,`],font=[TIMES,BOLD,14] ):
txtpp := textplot([0.24, 11, `l = 5`], font=[SYMBOL, 14]):
pntP := plot([[0.29, 10.87],[0.38, evalf(subs(t=0.38, Mp))]],color=black):
Pplot := display({Pp, txtp, txtpp, pntP}):
n := 10:
p := 1/2: # p = lambda/n
Pb10 := plot(Mb, t = -.4 .. .4, xtickmarks=xtics, ytickmarks=ytics, color=blue, thickness=2):
txtb10 := textplot([0.12, 6.4, `b(10, 1/2)`], font=[TIMES,BOLD,14]):
pntb10 := plot([[0.21, 6.3],[0.34, evalf(subs(t=0.34, Mb))]],color=black):
B[1] := display({Pb10, txtb10, pntb10, Pplot}):
n := 20:
p := 1/4:
Pb20 := plot(Mb, t = -.4 .. .4, xtickmarks=xtics, ytickmarks=ytics, color=blue, thickness=2):
txtb20 := textplot([0.12, 8, `b(20, 1/4)`], font=[TIMES,BOLD,14]):
pntb20 := plot([[0.214, 7.88],[0.35, evalf(subs(t=0.35, Mb))]],color=black):
B[2] := display({Pb20, txtb20, pntb20, Pplot}):
n := 50:
p := 1/10:
Pb50 := plot(Mb, t = -.4 .. .4, xtickmarks=xtics, ytickmarks=ytics, color=blue, thickness=2):
txtb50 := textplot([0.122, 9.5, `b(50, 1/10)`], font=[TIMES,BOLD,14]):
pntb50 := plot([[0.225, 9.4],[0.36, evalf(subs(t=0.36, Mb))]],color=black):
B[3] := display({Pb50, txtb50, pntb50, Pplot}):
display([seq(B[k], k = 1 .. 3)], insequence = true, title=`Convergence of moment-generating functions`,titlefont=[TIMES,BOLD,14]);

[Maple Plot]

> display([seq(B[k], k = 1 .. 3)], title=`Convergence of moment-generating functions`,titlefont=[TIMES,BOLD,14]);

[Maple Plot]

Now do the same with the respective probability histograms.

> xtics := [k$k=1..15]:

> P1 := ProbHistB(PoissonPDF(5, x), 0 .. 15):
txtp := textplot([11.2,0.062, `Poisson,`], font=[TIMES, BOLD, 14], color=black):
txtpl := textplot([13.9, 0.062, `l = 5`], font=[SYMBOL, 14]):
txtpntp := plot([[9.5,0.06], [8.5, 0.06]], color=black):
P2 := ProbHistFill(BinomialPDF(10, 1/2, x), 0 .. 15):
txtb := textplot([11.5, 0.18, `b(10, 1/2)`], font=[TIMES,BOLD,14]):
txtpntb := plot([[9.2, 0.177],[6.5, 0.177]], color=black):
P3 := plot([[0,0],[0,0]], x = -0.5 .. 15.9, xtickmarks=xtics):
H[1] := display({P1, txtp, txtpl, txtpntp, P2, txtb, txtpntb,P3}):
txtp := textplot([11.5,0.065, `Poisson,`], font=[TIMES, BOLD, 14], color=black):
txtpl := textplot([14.3, 0.064, `l = 5`], font=[SYMBOL, 14]):
txtpntp := plot([[9.6,0.064], [8.5, 0.064]], color=black):
P2 := ProbHistFill(BinomialPDF(20, 1/4, x), 0 .. 15):
txtb := textplot([11.5, 0.168, `b(20, 1/4)`], font=[TIMES,BOLD,14]):
txtpntb := plot([[9.4, 0.165],[6.5, 0.165]], color=black):
H[2] := display({P1,txtp, txtpl, txtpntp, P2,txtb, txtpntb,P3}):
txtp := textplot([12.2,0.037, `Poisson,`], font=[TIMES, BOLD, 14], color=black):
txtpl := textplot([14.9, 0.037, `l = 5`], font=[SYMBOL, 14]):
txtpntp := plot([[10.5,0.036], [9.5, 0.036]], color=black):
P2 := ProbHistFill(BinomialPDF(50, 1/10, x), 0 .. 15):
txtb := textplot([11.5, 0.153, `b(50, 1/10)`], font=[TIMES,BOLD,14]):
txtpntb := plot([[9.2, 0.152],[6.5, 0.152]], color=black):
H[3] := display({P1,txtp, txtpl, txtpntp, P2,txtb, txtpntb,P3}):
display([seq(H[k], k = 1 .. 3)], insequence = true, title=`Comparison of Probability Histograms`, titlefont=[TIMES,BOLD,14]);

[Maple Plot]

> f := PoissonPDF(lambda,x);

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

> lambda := 'lambda':
eta := lambda -> sum(t^x*f, x = 0 .. infinity);

eta := proc (lambda) options operator, arrow; sum(t...

> eta(lambda);
series(eta(lambda), t, 16);

exp(-5)*exp(5*t)

series(exp(-5)+5*exp(-5)*t+25/2*exp(-5)*t^2+125/6*e...
series(exp(-5)+5*exp(-5)*t+25/2*exp(-5)*t^2+125/6*e...
series(exp(-5)+5*exp(-5)*t+25/2*exp(-5)*t^2+125/6*e...

> lambda := 5;
evalf(series(eta(lambda), t, 16));

lambda := 5

series(.6737946999e-2+.3368973500e-1*t+.8422433749e...
series(.6737946999e-2+.3368973500e-1*t+.8422433749e...
series(.6737946999e-2+.3368973500e-1*t+.8422433749e...

> eta[b] := sum(t^x*BinomialPDF(n, p, x), x = 0 .. n):

> n := 50;
p := 1/10;
evalf(series(eta[b], t, 16));

n := 50

p := 1/10

series(.5153775207e-2+.2863208449e-1*t+.7794289665e...
series(.5153775207e-2+.2863208449e-1*t+.7794289665e...
series(.5153775207e-2+.2863208449e-1*t+.7794289665e...

Return to Menu

2. Convergence of moment-generating functions to illustrate the Central Limit Theorem -- Exponential Distribution

Take a random sample of size n from an exponential distribution with theta = 2. Find the moment-generating function of W = (Xbar-theta)/(theta/sqrt(n)) and show that it converges to that for the N (0, 1) distribution.

> n := 'n':
assume(t < 1/2):
ME := simplify(int(exp(x*t)*ExponentialPDF(2,x), x=0..infinity)); #Exp mgf

ME := -1/(2*t-1)

> M := (subs(t = t/2/sqrt(n), ME))^n*exp(-sqrt(n)*t); # W mgf

M := (-1/(t/(sqrt(n))-1))^n*exp(-sqrt(n)*t)

> #M := exp(-sqrt(n)*t)/(1 - t/sqrt(n))^n;
MZ := int(exp(z*t)*NormalPDF(0, 1, z), z = -infinity .. infinity);

MZ := exp(1/2*t^2)

> xtics := [-1+0.2*k$k=0..10]:
ytics := [0.8+0.1*k$k=0..10]:

> for m from 1 to 15 do
n := 15*m:
PZ := plot(MZ, t = -1 .. 1, y=0.71.. 1.82, xtickmarks=xtics, ytickmarks=ytics, color=black, thickness = 2, labels=[``,``]):
PU := plot(M, t = -1 .. 1, y = 0.71 .. 1.82, xtickmarks=xtics, ytickmarks=ytics, color=blue, thickness=2, labels=[``,``]):
txt := textplot([0.4,1.5,`n = `||n], font=[TIMES,BOLD,14]):
H[m] := display({PZ, PU, txt}):
od:
display([seq(H[m], m = 1 .. 15)], insequence=true, title=`Limit of Exponential Moment-generating\nFunctions to that of the N(0,1) [in black]`, titlefont=[TIMES,BOLD,14]);

[Maple Plot]



Sum of chi^2 (1) Random Variables

Let X[1] , X[2] , ..., X[n] be a random sample of size n from a chi^2 (1) distribution. Then the distribution of Y = sum(X[i],i = 1 .. n) is chi^2 ( n ). If we let W = (Y-n)/sqrt(2*n) , then E( W ) = 0 and Var( W ) = 1. The following animation shows that the p.d.f. of W approaches the p.d.f. for the N (0, 1) distribution as n increases. However, because of the skewness of the chi^2 (1) distribution, n must be quite large for W to have an approximate N (0, 1) distribution.

> k := 'k':
xtics := [-3.5+0.5*k$k=0..14]:
ytics := [0.05*k$k=1..8]:

> w := 'w':
n := 'n':
normplot := plot(NormalPDF(0, 1, x), x = -3.52 .. 3.52, y = 0 .. 0.43, xtickmarks=xtics, ytickmarks=ytics, color = blue, thickness=2):
for k from 1 to 15 do
n := 20*k:
f := ChisquarePDF(n, x):
h[k] := changevar(x = sqrt(2*n)*w + n, sqrt(2*n)*f):
P[k] := plot(h[k], w = -3.52 .. 3.52, xtickmarks=xtics, ytickmarks=ytics, color=red, thickness=2, labels=[``,``]):
txt := textplot([1.8, 0.3, `n = `||n], font=[TIMES,BOLD,14]):
H[k] := display({normplot, P[k], txt}):
od:

> display([seq(H[k], k = 1 .. 15)], title = `p.d.f. of the Sum of Transformed Chi-square(1)\nRandom Variables and the N(0,1) p.d.f.`, insequence = true, titlefont=[TIMES,BOLD,14]);

[Maple Plot]

Return to Menu

3. Convergence of moment-generating functions to illustrate the Central Limit Theorem -- A U-shaped Distribution

> g := NormalPDF(0, 1, x);
MZ := int(exp(t*x)*g, x = -infinity .. infinity);

g := 1/2*sqrt(2)*exp(-1/2*x^2)/(sqrt(Pi))

MZ := exp(1/2*t^2)

> f := 3/2*x^2;

f := 3/2*x^2

> MU := int(exp(t*x)*f, x = -1 .. 1);

MU := 3/2*(exp(t)*t^2+2*exp(t)-2*exp(t)*t-exp(-t)*t...

> limit(MU,t=0,left);
limit(MU,t=0,right);

1

1

> MUp := diff(MU, t);

MUp := 3/2*(exp(t)*t^2+exp(-t)*t^2)/(t^3)-9/2*(exp(...

> mu := limit(MUp, t=0, left);
mu := limit(MUp, t=0, right);

mu := 0

mu := 0

> MUpp := diff(MUp, t);

MUpp := 3/2*(exp(t)*t^2+2*exp(t)*t-exp(-t)*t^2+2*ex...
MUpp := 3/2*(exp(t)*t^2+2*exp(t)*t-exp(-t)*t^2+2*ex...

> var := limit(MUpp, t=0, left);
var := limit(MUpp, t=0, right);

var := 3/5

var := 3/5

> n := 'n':
MUn := (subs(t = sqrt(5/n/3)*t, MU))^n:

> MUnp := diff(MUn,t):

> limit(MUnp, t=0, left);
limit(MUnp, t=0, right);

0

0

> MUnpp := diff(MUnp, t):

> limit(MUnpp, t=0, left);
limit(MUnpp, t=0, right);

1

1

> n := 'n':
MUn := (subs(t = sqrt(5/n/3)*t, MU))^n;

MUn := (9/50*(5/3*exp(1/3*sqrt(15)*sqrt(1/n)*t)*t^2...
MUn := (9/50*(5/3*exp(1/3*sqrt(15)*sqrt(1/n)*t)*t^2...
MUn := (9/50*(5/3*exp(1/3*sqrt(15)*sqrt(1/n)*t)*t^2...

> k := 'k':
xtics := [-2 + 0.5*k$k=0 .. 8]:
ytics := [0.5*k$k=1..7]:

> for n from 1 to 12 do
PZ := plot(MZ, t = -2 .. 2, y=0.45.. 3.55, xtickmarks=xtics, ytickmarks=ytics, color=black, thickness=2, labels=[``,``]):
PU := plot(MUn, t = -2 .. 2, y = 0.5 .. 3.5, xtickmarks=xtics, ytickmarks=ytics, color=blue, thickness=2, labels=[``,``]):
txt := textplot([0.5,2.5,`n = `||n],font=[TIMES,BOLD,14]):
H||n := display({PZ, PU, txt}):
od:
display([seq(H||n, n = 1 .. 12)], insequence=true, title=`Graphs of Moment-Generating Functions for\nX-bar from U-shaped Distribution and N(0,1)`, titlefont=[TIMES,BOLD,14]);

[Maple Plot]

Empirical Support for Convergence

Theoretical Support for Convergence

Return to Menu