Elliot Tanis
June 12, 2006
Figures 5.4-2(a) and 5.4-2(b)  (Fig5_4-2.mws)

 > randomize():

 > for k from 1 to 1000 do X := [seq(-ln(1 - rng()), k = 1 .. 16)]: TT[k] := (Mean(X) - 1)/(StDev(X)/4): od: T := [seq(TT[k], k = 1 .. 1000)]:

 > tbar := Mean(T);

 > s := StDev(T);

 > Min(T), Max(T);

 > xtics := [seq(-7 + k*0.5, k = 0 .. 21)]: ytics := [seq(0.05*k, k = 1 .. 8)]: P1 := HistogramFill(T, -6.5 .. 3.5, 20): P2 := plot([[0,0],[0,0]], x = -6.5 .. 3.5, y = 0 .. 0.41, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): display({P1, P2});

Similar to Figure 5.4-2(a)

 > #zper := [seq(NormalP(0, 1, k/1001), k = 1 .. 1000)]: k := 'k': for k from 1 to 500 do   zp[k] := NormalP(0, 1, k/1001):   zp[1000-k+1] := -zp[k]: od:

 > zper := [seq(zp[k], k = 1 .. 1000)]:

 > Torder := sort(T):

 > ScatPlotPoint(Torder, zper);

Similar to Figure 5.4-2(b)

 >

 > X := [seq(-ln(1 - rng()) + 9, k = 1 .. 16)]; Mean(X); StDev(X);

The above are newly generated data. The following are those data in the text. We use the following at this time.

 > X := [11.9776, 9.3889, 9.9798, 13.4676, 9.2895, 10.1242, 9.5798, 9.3148, 9.0605, 9.1680, 11.0394, 9.1083, 10.3720, 9.0523, 13.2969, 10.5852];

 > xbar := Mean(X); sx := StDev(X);

 > probs := [seq(1/16, k = 1 .. 16)]; empPDF := zip((x,y)->(x,y),X,probs);

 > for k from 1 to 1000 do XX := DiscreteS(empPDF, 16): TT[k] := (Mean(XX) - 10.3003)/(StDev(XX)/4): od: Tre := [seq(TT[k], k = 1 .. 1000)]: tbar := Mean(Tre); s := StDev(Tre); Min(Tre), Max(Tre);

 > xtics := [seq(-6.5+0.5*k, k = 0 .. 20)]: ytics := [seq(0.05*k, k = 1 .. 8)]:

 > P1 := HistogramFill(Tre, -6.5 .. 3.5, 20): P2 := plot([[0,0],[0,0]], x=-6.5 .. 3.5, y = 0 .. 0.445, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): display({P1, P2});

Similar to Figure 5.4-3(a)

 > Treorder := sort(Tre):

 > ScatPlotPoint(Treorder, zper);

Similar to Figure 5.4-3(b)

Now find the 0.025 and 0.975 quantiles

 > c := Percentile(Tre, 0.025); d := Percentile(Tre, 0.975);

 > left := evalf(xbar - d*sx/4);

 > right := evalf(xbar - c*sx/4);

 > xbar;

 > sx;

Hopefully "10" is in the interval

 > [left, right];

 >