•  1.    Simple Simulations
•  2.    Relationships Between Distributions
•  3.    The Power Function of a Statistical Test
•  4.    Minimum-Length Confidence Intervals for Variances
•  5.    Investigating the Properties of a Distribution
•  6.    Understanding an Important Theorem
•  7.    Relationships Between Random Variables
•  8.    Confidence Intervals
•  9.    Order Statistics
•  10.  Regression

Examples for Statistics Supplement

by Zaven Karian

 1.    Simple Simulations

>	restart: with(plots, display): libname:="C:/mylib/stat",libname: with(stat):

>	A := Die(6,8);

>	Freq(A, 1..6);

>	A := Die(4, 400):

>	Freq(A, 1..4);

>	Coins := DiscreteS( [0, 1/2, 1, 1/2], 8);

>	pdf := x/10;

>	Sample := DiscreteS(pdf, 1..4, 5);

>	Sample := DiscreteS(pdf, 1..4, 300):

>	EH := Histogram(Sample, 0.5..4.5, 4):

>	PH := ProbHist(pdf, 1..4):

>	display({EH,PH});

>	S1 := BinomialS(10, 0.25, 200):

>	S2 := GammaS(3, 2, 100):

>	PlotEmpPDF(S2);

>	PlotEmpPDF(S2, 0..25, 10);

>	PlotEmpCDF(S2,0..25);

>	BoxWhisker(S1,S2);

>	QQ(S1,S2);

>

 2.    Relationships Between Distributions

>	restart: with(plots,display): libname:="C:/mylib/stat",libname: with(stat):

>	for i from 1 to 16 do

>	H := ProbHist(BinomialPDF(2*i, 0.25, x), 0..16, 17):

>	N := plot(NormalPDF(i/2, 3*i/8, x), x=i/2-4*sqrt(3*i/8)..i/2+4*sqrt(3*i/8)):

>	P||i := display({H,N}):

>

od:

>	display([seq(P||i, i=1..16)], insequence = true);

>	display(P2);

>	display(P3);

>	display(P12);

>	display(P16);

>

 3.    The Power Function of a Statistical Test

>	restart: with(plots,display): libname:="C:/mylib/stat",libname: with(stat):

>	n := 21; var := 100; alpha := 0.05;

>	L := var*ChisquareP(n-1, alpha/2)/sigma^2;

>	U := var*ChisquareP(n-1, 1-alpha/2)/sigma^2;

>	K := 1-ChisquareCDF(n-1,U) + ChisquareCDF(n-1,L);

>	P4 := plot(K, sigma=0..20): P4;

>

>	display(P4);

>	Alpha := [0.005, 0.01, 0.025, 0.05, 0.1, 0.2];

>	for i from 1 to 6 do

>	alpha := Alpha[i];

>	L := var*ChisquareP(n-1, alpha/2)/sigma^2;

>	U := var*ChisquareP(n-1, 1-alpha/2)/sigma^2;

>	K := 1-ChisquareCDF(n-1,U) + ChisquareCDF(n-1,L):

>	P||i := plot(K, sigma=0..20):

>

od:

>	display({P1,P2,P3,P4,P5,P6});

>	display([P1,P2,P3,P4,P5,P6], insequence = true);

>

>	alpha := 0.05;

>	for i from 1 to 6 do

>	n := 1+4*i:

>	L := var*ChisquareP(n-1,alpha/2)/sigma^2:

>	U := var*ChisquareP(n-1,1-alpha/2)/sigma^2:

>	K := 1-ChisquareCDF(n-1,U) + ChisquareCDF(n-1,L):

>	Q||i := plot(K, sigma=0..20):

>

od:

>	display({Q1,Q2,Q3,Q4,Q5,Q6});

>	display([Q1,Q2,Q3,Q4,Q5,Q6], insequence = true);

>

 4.    Minimum-Length Confidence Intervals for Variances

Minimal length CI for variance.

Using the ordianry CI (with equal probability in each tail),

we get the following for the length of the CI for the variance.

 

>	restart: with(plots, display):  libname:="C:/mylib/stat",libname: with(stat):

>	n := 7: alpha := 0.05:

>	UsualL := S^2*(n-1)*( 1/ChisquareP(n-1,alpha/2) - 1/ChisquareP(n-1,1-alpha/2) );

>	f:=ChisquarePDF(n-1,x);

>	eq1 := int(f, x=a..b)=1-alpha;

>	eq2:=a^2*subs(x=a,f)=b^2*subs(x=b,f);

>	solution := fsolve({eq1,eq2},{a,b},{a=0..n-1, b=n-1..infinity});

>	assign(%);

>	MinimumL := S^2*(n-1)*(1/a-1/b);

>	100*(UsualL-MinimumL)/UsualL;

>

>	restart: with(plots, display): libname:="C:/mylib/stat",libname: with(stat):

>

n := 7:

>	UsualL:=(n-1)*(1/ChisquareP(n-1,alpha/2)-1/ChisquareP(n-1,1-alpha/2));

>	A := plot(UsualL, alpha=0.01..0.2):

>	f:=ChisquarePDF(n-1,x):

>	eq1 := int(f, x=a..b)=1-alpha:

>	eq2:=a^2*subs(x=a,f)=b^2*subs(x=b,f):

>	MinLenArray := array(1..39):

>	for i from 2 to 40 do

>	eqns := subs(alpha=0.005*i, {eq1, eq2}); Sol := fsolve(eqns, {a,b}, {a=0..n-1, b=n-1..infinity});

>	assign(Sol);

>	MinLenArray[i-1] := [0.005*i, (n-1)*(1/a-1/b)];

>	a:='a'; b:='b';

>

od:

>	MinLenList := convert(MinLenArray, list);

>	A := plot(UsualL, alpha=0.01..0.2):

>	B:=plot(MinLenList, color=blue):

>	display({A,B});

>	Improv := [seq([0.005+0.005*i,(value(subs(alpha=0.005+0.005*i,UsualL)-MinLenList[i][2]))], i=1..39)];

>	PercImprov := [seq([0.005+0.005*i, 100*(value(subs(alpha=0.005+0.005*i,UsualL))-MinLenList[i][2])/value(subs(alpha=0.005+0.005*i,UsualL))], i=1..39)];

>	StandLen := S^2*(n-1)*( 1/ChisquareP(n-1,alpha/2) - 1/ChisquareP(n-1,1-alpha/2) );

>	LFact := subs({alpha=0.05, S=1}, StandLen);

>	PlotList1 := [seq( [i, value(subs(n=i,LFact))], i=3..15)];

>	A := plot(PlotList1, color=red): A;

>	f := ChisquarePDF(n-1, x);

>	eq1 := int(f, x=a..b)=1-alpha;

>	eq2 := a^2*subs(x=a,f) = b^2*subs(x=b,f);

>	alpha := 0.05;

>	PlotArray := array(1..13);

>	for i from 3 to 15 do

>	eqns := subs(n=i, {eq1, eq2}); Sol := fsolve(eqns, {a,b}, {a=0..i-1, b=i-1..infinity});

>	assign(Sol);

>	PlotArray[i-2] := [i, (i-1)*(1/a-1/b)];

>	a:='a'; b:='b';

>

od:

>	PlotList2 := convert(PlotArray, list);

>	B:=plot(PlotList2, color=blue):

>	display({A,B});

>	[seq( [PlotList1[i][1], 100*(PlotList1[i][2] - PlotList2[i][2])/PlotList1[i][2]], i=1..13)];

>

 5.    Investigating the Properties of a Distribution

>	restart:  with(plots, display, animate): libname:="C:/mylib/stat",libname: with(stat):

>	f := ChisquarePDF(nu, x);

>	assume(nu>0); interface(showassumed=0);

>	simplify(int(f, x=0..infinity));

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

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

>	var := simplify(%);

>	animate(f, x=0..25, nu=1..15);

>	plot3d(f, x=0..10, nu=0..5, orientation=[135,45], color=blue, axes=boxed);

>	df := simplify(diff(f, x));

>	solve(df=0, x);

>

 6.    Understanding an Important Theorem

>	restart: with(plots, display):  libname:="C:/mylib/stat",libname: with(stat):

>	f := x -> (3/2)*x^2;

>	int(f(x), x=-1..1);

>	mu := int(x*f(x), x=-1..1);

>	var := int(x^2*f(x), x=-1..1)- mu^2;

>	F := int(f(t), t=-1..x);

>	A := ContinuousS(f(x), -1..1, 5);

>	A := ContinuousS(f(x), -1..1, 400):

>	H := Histogram(A, -1..1, 12):

>	P := plot(f(x), x=-1..1):

>	display({P,H});

>	S4 := [seq(ContinuousS(f(x), -1..1, 4), i=1..300)]:

>

S4[1];

>	M4 := [seq(Mean(S4[i]), i=1..300)]:

>	H4 := Histogram(M4, -1..1, 7):

>	n4 := NormalPDF(mu, var/4, x):

>	N4 := plot(n4, x=-1..1, color=blue):

>	display({P, H4, N4});

>	S8 := [seq(ContinuousS(f(x), -1..1, 8), i=1..300)]:

>

A[1];

>	M8 := [seq(Mean(S8[i]), i=1..300)]:

>	H8 := Histogram(M8, -1..1, 7):

>	n8 := NormalPDF(mu, var/8, x):

>	N8 := plot(n8, x=-1..1, color=blue):

>	display({P, H8, N8});

>

>	A := [seq(ContinuousS(f(x), -1..1, 16), i=1..300)]:

>	M16 := [seq(Mean(A[i]), i=1..300)]:

>	H16 := Histogram(M16, -1..1, 13):

>	n16 := NormalPDF(mu, var/16, x):

>	N16 := plot(n16, x=-1..1, color=blue):

>	display({P, H16, N16});

>

 7.    Relationships Between Random Variables

>	restart: with(plots, display): libname:="C:/mylib/stat",libname: with(stat):

>	Z1 := NormalS(0,1,500):

>	Z1H := Histogram(Z1, -4..4, 9):

>	N01 := plot(NormalPDF(0,1,x), x=-4..4):

>	display({Z1H, N01});

>	Y := [seq(Z1[i]^2, i=1..500)]:

>	YH := Histogram(Y, 0..12, 16):

>	CH1 := plot(ChisquarePDF(1,x), x=0..12):

>	display({YH, CH1});

>

>	Z2 := NormalS(0,1,500):

>	Z3 := NormalS(0,1,500):

>	Y := [seq(Z1[i]^2+Z2[i]^2+Z3[i]^2, i=1..500)]:

>	YH := Histogram(Y, 0..18, 24):

>	CH3 := plot(ChisquarePDF(3,x), x=0..12):

>	#interface(plotdevice=postscript, plotoutput=Fig7b);

>	display({YH, CH3});

>

>	Z := [seq(NormalS(37, 16, 3), i=1..200)]:

>	S := [seq(sum( (Z[i][j]-37)^2/16, j=1..3), i=1..200)]:

>	T := [seq(sum( (Z[i][j]-Mean(Z[i]))^2/16, j=1..3), i=1..200)]:

>	SH := Histogram(S, 0..20, 16):

>	TH := Histogram(T, 0..20, 16):

>	CH2 := plot(ChisquarePDF(2,x), x=0..20):

>	CH3 := plot(ChisquarePDF(3,x), x=0..20):

>	display({CH3, SH});

>	display({CH2, TH});

>

 8.    Confidence Intervals

>	restart: with(plots,display): libname:="C:/mylib/stat",libname: with(stat):

>	ListOfSamples := [seq(NormalS(40,12,5), i=1..50)]:

>	CIVarUnknown := ConfIntMean(ListOfSamples, 80):

>	ConfIntPlot(CIVarUnknown, 40);

>	CIVarKnown := ConfIntMean(ListOfSamples, 80, 12):

>	ConfIntPlot(CIVarKnown, 40);

>

 9.    Order Statistics

>	restart: with(plots, display): libname:="C:/mylib/stat",libname: with(stat):

A study of the p.d.f.s of order stat

First, take example 10.1-3

>	f := x -> 1;

>	F := x -> int(f(t), t=0..x);

>	g := (n, r, y) -> (n!/((r-1)!*(n-r)!)) * F(y)^(r-1)*(1-F(y))^(n-r)*f(y);

>	G:= (n,r,y) -> sum(n!/(k!*(n-k)!)*F(y)^k*(1-F(y))^(n-k), k=r..n);

>	pdf := g(4,3,y);

>	int(pdf, y=1/3..2/3);

>	G(4,3,2/3)-G(4,3,1/3);

Note that this "hides the expressions for g and G.

We get around this by the following.

Can this be done if we assume that the underlying distribution

is more complicated? Assume we have an exponential.

>	assume(theta > 0); interface(showassumed=0);

>	f := x-> ExponentialPDF(theta, x);

>	F := x-> int(f(t), t=0..x);

>	g := (n, r, y) -> (n!/((r-1)!*(n-r)!)) * F(y)^(r-1)*(1-F(y))^(n-r)*f(y);

>	G:= (n,r,y) -> sum(n!/(k!*(n-k)!)*F(y)^k*(1-F(y))^(n-k), k=r..n);

We can now compute probabilities. For example, if n=4,

what is the probability that the third order statistic is

between 1/2 and 1.

>	P := simplify(G(4,3,1)-G(4,3,1/2));

If we assume a specific theta such as theta=1, we get

a numeric answer.

>	evalf(subs(theta=1,P));

We can now study the effect of theta on an order

statistic for specific values of n and r. Let's

take n=4, r=3.

>	pdf := simplify(g(4,3,y));

>

N := 5;

>	for i from 1 to N do

>	gg := simplify(subs(theta=i, g(4,3,y)));

>	P||i := plot(gg, y=0..10, color=blue):

>

od:

>	display([seq(P||i, i=1..5)], insequence=true);

For a fixed theta, say theta=1, we can also consider

the shapes of the order stat for r=1, 2, 3, ...

Taking theta =1 and n=4, we look at r=1,2,3, and 4.

>	for i from 1 to 4 do

>	gg := simplify(subs(theta=1, g(4,i,y)));

>	P||i := plot(gg, y=0..4, color=blue):

>

od:

>	display({P1,P2,P3,P4});

Unfortunately, the expressions for g and G in both examples

"hide" their actual forms. The following will produce more

explicit representations of g and G, if that is desired.  

>	assume( theta > 0); interface(showassumed=0);

>	f := ExponentialPDF(theta, y);

>	int(f, y=0..t);

>	F := subs(t=y,%);

>	g := (n!/((r-1)!*(n-r)!)) * F^(r-1)*(1-F)^(n-r)*f;

>	G := sum(n!/(k!*(n-k)!)*F^k*(1-F)^(n-k), k=r..n);

For contrasrt, we now consider order stat from a

symmetric, two-parameter distribution: N(mu,var).

>	assume(sigma > 0); interface(showassumed=0);

>	f := NormalPDF(mu, sigma^2, y);

>	int(f, y=-infinity..t);

>	F := subs(t=y,%);

>	g := (n!/((r-1)!*(n-r)!)) * F^(r-1)*(1-F)^(n-r)*f;

>	G := simplify(sum(n!/(k!*(n-k)!)*F^k*(1-F)^(n-k), k=r..n));

As before, we can calculate probabilities in

specific cases.

>	GG := simplify(subs({n=4, r=3, mu=0, sigma=1}, G));

>	evalf(subs(y=1/2, GG)- subs(y=-1/2, GG));

And we can obtain expressions and graphs for

specific p.d.f.s of order stat or through

animation observe the effect of the symmetry

of the underlying distribution by comparing

g_i with g_(n-i+1).

>	gg := simplify(subs({n=4, r=3, mu=0, sigma=1}, g));

>	plot(gg, y=-3.5..3.5);

>

N := 10;

>	for i from 1 to N do

>	gg := simplify(subs({n=N, r=i, mu=0, sigma=1},g));

>	P||i := plot(gg, y=-3.5..3.5, color=blue):

>

od:

>	display([seq(P||i, i=1..N)], insequence=true);

>

N:=20;

>	for i from 1 to N/2 do

>	gg := simplify(subs({n=N, r=i, mu=0, sigma=1},g));

>	ggg := simplify(subs({n=N, r=N-i+1, mu=0, sigma=1},g)); P||i := plot({gg,ggg}, y=-3.5..3.5, color=blue):

>

od:

>	display([seq(P||i, i=1..N/2)], insequence=true);

>

 10.  Regression

>	restart: with(plots, display):  libname:="C:/mylib/stat",libname: with(stat):

>	X := [1,2,3,4,5];

>	Y := [2,3,1,5,4];

>	r:= Correlation(X,Y);

>	A := BivariateNormalS(100, 225, 400, 625, 0.9, 200):

>

A[1..5];

>	r := Correlation(A);

>	y := LinReg(A, x);

>	Y := PolyReg(A, 3, x);

>	ScatPlot(A);

>	display( ScatPlot(A), plot(y, x=50..150));

>