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

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

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

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

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

> f := PoissonPDF(lambda,x);

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

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

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

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

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

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 Find the moment-generating function of 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

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

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

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

Let , , ..., be a random sample of size n from a (1) distribution. Then the distribution of is ( n ). If we let , 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 (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]);

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

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

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

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

> MUp := diff(MU, t);

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

> MUpp := diff(MUp, t);

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

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

> MUnpp := diff(MUnp, t):

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

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

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