Example 4: Minimum-Length Confidence Intervals for Variances

by Zavan Karian

> restart: with(plots, display): libname:="C:/mylib/statistics",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;

> n := 7:

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

Warning, incomplete string;  use " to end the string

> plot(UsualL, 0.01..0.2, 4.4..4.6);

> 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, 0.01..0.2):

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

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

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