Elliot Tanis
April 20, 2006
Exercise 5.3-12  (Exercise_5_3-12.mws)

(a) Simulations of T  when sampling from N (1, 1).

 > randomize(): for k from 1 to 500 do   X := NormalS(1,1,7):   TT[k] := evalf((Mean(X) - 1)/(StDev(X)/sqrt(7))): od: T := [seq(TT[k], k = 1 .. 500)]: Min(T),Max(T); tbar := Mean(T); mu := 0; svar := Variance(T); var := 6/4;

 > xtics := [seq(-5 + 0.5*k, k = 0 .. 20)]: ytics := [seq(0.05*k, k = 1 .. 9)]: P1 := plot(TPDF(6,x), x = -5 .. 5, color=black, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): P2 := HistogramFill(T, -5.5 .. 5.5, 22): display({P1, P2});

 > quantiles := [seq(TP(6, k/501), k = 1 .. 500)]:

 > T := sort(T):

 > xtics := [seq(-5.5 + k*.5, k = 0 .. 22)]: ytics := xtics: P1 := plot([[-5,-5],[5,5]], x = -5.5 .. 5.5, y = -5.5 .. 5.5, color=black, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): P2 := ScatPlotCirc(T, quantiles): display({P1, P2});

(b) Now sample from an exponential distribution with mean 1 and observe values of " T ".

 > for k from 1 to 500 do   X := ExponentialS(1,7):   TT[k] := evalf((Mean(X) - 1)/(StDev(X)/sqrt(7))): od: T := [seq(TT[k], k = 1 .. 500)]: Min(T),Max(T); tbar := Mean(T); mu := 0; svar := Variance(T); var := 6/4;

 > xtics := [seq(-7 + 0.5*k, k = 0 .. 28)]: ytics := [seq(0.05*k, k = 1 .. 9)]: P1 := plot(TPDF(6,x), x = -7 .. 7, color=black, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): P2 := HistogramFill(T, -7 .. 7, 28): display({P1, P2});

If you have already done the previous simulation, you do not have to again calculate the quantiles for the t-distribution with with 6 degrees of freedom. If this is the first simulation, remove # and run the following line.

 > #quantiles := [seq(TP(6, k/501), k = 1 .. 500)]:

 > T := sort(T):

 > xtics := [seq(-10 + k, k = 0 .. 18)]: ytics := xtics: P1 := plot([[-7,-7],[7,7]], x = -7 .. 7, y = -7 .. 7,color=black, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): P2 := ScatPlotCirc(T, quantiles): display({P1, P2});

(c) Now try sampling from the Cauchy distribution.  The following procedure can be used.

 > CauchyS(7);

 > for k from 1 to 500 do   X := CauchyS(7):   TT[k] := evalf(Mean(X)/(StDev(X)/sqrt(7))): od: T := [seq(TT[k], k = 1 .. 500)]: Min(T),Max(T); tbar := Mean(T); mu := 0; svar := Variance(T); var := 6/4;

 > xtics := [seq(-7 + 0.5*k, k = 0 .. 28)]: ytics := [seq(0.05*k, k = 1 .. 9)]: P1 := plot(TPDF(6,x), x = -4 .. 4, color=black, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): P2 := HistogramFill(T, -4 .. 4, 16): display({P1, P2});

 > T := sort(T):

 > xtics := [seq(-5 + k*.5, k = 0 .. 20)]: ytics := xtics: P1 := plot([[-4,-4],[4,4]], x = -4 .. 4, y = -4 .. 4, color=black, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): P2 := ScatPlotCirc(T, quantiles): display({P1, P2});

 >