>    restart:
read `C:\\Tanis-Hogg\\Maple Examples\\stat.m`:
read `C:\\Tanis-Hogg\\Maple Examples\\HistogramFill.txt`:
read `C:\\Tanis-Hogg\\Maple Examples\\ScatPlotCirc.txt`:
read `C:\\Tanis-Hogg\\Maple Examples\\ScatPlotPoint.txt`:
with(plots):

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

Figure 5.4-1(a) Histogram of the Trimmed Means and Figure 5.4-1(b) N(0, 1) Quantiles versus Trimmed Means Quantiles

Using the values of W  listed in Section 5.4, the following N  = 1000 values of the trimmed mean were simulated. That is, 1000 times we took a resample and calculated the value of the trimmed mean, the mean of the middle 10 order statistics, yielding the following values.

>    X := [4.30, 4.33, 4.35, 4.40, 4.44, 4.46, 4.46, 4.47, 4.47, 4.48, 4.49, 4.49, 4.50, 4.50, 4.51,4.53, 4.54, 4.54, 4.55, 4.55, 4.57, 4.57, 4.58, 4.58, 4.58, 4.58, 4.58, 4.59, 4.60, 4.60,4.60, 4.61, 4.61, 4.61, 4.61, 4.61, 4.61, 4.62, 4.62, 4.62, 4.63, 4.63, 4.63, 4.64, 4.64,4.64, 4.64, 4.64, 4.64, 4.64, 4.65, 4.65, 4.66, 4.67, 4.67, 4.67, 4.67, 4.67, 4.67, 4.67,4.67, 4.67, 4.67, 4.67, 4.67, 4.67, 4.67, 4.67, 4.67, 4.68, 4.68, 4.68, 4.69, 4.69, 4.69,4.69, 4.69, 4.69, 4.69, 4.70, 4.70, 4.70, 4.70, 4.70, 4.70, 4.70, 4.71, 4.71, 4.71, 4.71,4.71, 4.71, 4.71, 4.71, 4.71, 4.71, 4.72, 4.72, 4.72, 4.72, 4.72, 4.72, 4.72, 4.72, 4.72,4.72, 4.72, 4.72, 4.72, 4.72, 4.73, 4.73, 4.73, 4.73, 4.73, 4.73, 4.73, 4.73, 4.73, 4.74,4.74, 4.74, 4.74, 4.74, 4.74, 4.74, 4.74, 4.74, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75, 4.75,4.75, 4.75, 4.76, 4.76, 4.76, 4.76, 4.76, 4.76, 4.76, 4.76, 4.77, 4.77, 4.77, 4.77, 4.77,4.77, 4.77, 4.77, 4.77, 4.77, 4.78, 4.78, 4.78, 4.78, 4.78, 4.78, 4.78, 4.78, 4.78, 4.79,4.79, 4.79, 4.79, 4.79, 4.79, 4.79, 4.79, 4.79, 4.79, 4.79, 4.79, 4.79, 4.80, 4.80, 4.80,4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.80, 4.81, 4.81,4.81, 4.81, 4.81, 4.81, 4.81, 4.81, 4.81, 4.81, 4.81, 4.81, 4.81, 4.81, 4.81, 4.82, 4.82,4.82, 4.82, 4.82, 4.82, 4.82, 4.82, 4.82, 4.82, 4.82, 4.82, 4.82, 4.82, 4.82, 4.82, 4.82,4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.84, 4.84, 4.84, 4.84,4.84, 4.84, 4.84, 4.84, 4.84, 4.84, 4.84, 4.84, 4.84, 4.84, 4.84, 4.84, 4.84, 4.84, 4.84,4.84, 4.84, 4.85, 4.85, 4.85, 4.85, 4.85, 4.85, 4.85, 4.85, 4.85, 4.85, 4.85, 4.85, 4.85,4.85, 4.86, 4.86, 4.86, 4.86, 4.86, 4.86, 4.86, 4.86, 4.86, 4.86, 4.86, 4.86, 4.86, 4.86,4.86, 4.86, 4.86, 4.86, 4.86, 4.86, 4.86, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87,4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87, 4.87,4.88, 4.88, 4.88, 4.88, 4.88, 4.88, 4.88, 4.88, 4.88, 4.88, 4.88, 4.88, 4.88, 4.88, 4.88,4.88, 4.88, 4.88, 4.88, 4.88, 4.88, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89,4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89,4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.90, 4.90, 4.90, 4.90,4.90, 4.90, 4.90, 4.90, 4.90, 4.90, 4.90, 4.90, 4.90, 4.90, 4.90, 4.90, 4.90, 4.90, 4.90,4.90, 4.90, 4.90, 4.91, 4.91, 4.91, 4.91, 4.91, 4.91, 4.91, 4.91, 4.91, 4.91, 4.91, 4.91,4.91, 4.91, 4.91, 4.91, 4.91, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92,4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.92, 4.93,4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93,4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 4.94, 4.94, 4.94,4.94, 4.94, 4.94, 4.94, 4.94, 4.94, 4.94, 4.94, 4.94, 4.94, 4.94, 4.94, 4.94, 4.94, 4.94,4.94, 4.94, 4.94, 4.94, 4.95, 4.95, 4.95, 4.95, 4.95, 4.95, 4.95, 4.95, 4.95, 4.95, 4.95,4.95, 4.95, 4.95, 4.95, 4.95, 4.95, 4.95, 4.95, 4.95, 4.95, 4.96, 4.96, 4.96, 4.96, 4.96,4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96,4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.96, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97,4.97, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97, 4.97,4.97, 4.97, 4.97, 4.97, 4.97, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98,4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98,4.98, 4.98, 4.98, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99,4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99,4.99, 4.99, 4.99, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00,5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.01,5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01,5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.01, 5.02, 5.02,5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02,5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.02, 5.03, 5.03, 5.03, 5.03,5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03,5.03, 5.03, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04,5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04, 5.04,5.04, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05,5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.05, 5.06, 5.06,5.06, 5.06, 5.06, 5.06, 5.06, 5.06, 5.06, 5.06, 5.06, 5.06, 5.06, 5.06, 5.06, 5.06,5.06, 5.06, 5.06, 5.07, 5.07, 5.07, 5.07, 5.07, 5.07, 5.07, 5.07, 5.07, 5.07, 5.07,5.07, 5.07, 5.07, 5.07, 5.07, 5.08, 5.08, 5.08, 5.08, 5.08, 5.08, 5.08, 5.08, 5.08,5.08, 5.08, 5.08, 5.08, 5.09, 5.09, 5.09, 5.09, 5.09, 5.09, 5.09, 5.09, 5.09, 5.09,5.09, 5.09, 5.09, 5.09, 5.09, 5.09, 5.10, 5.10, 5.10, 5.10, 5.10, 5.10, 5.10, 5.10,5.10, 5.10, 5.11, 5.11, 5.11, 5.11, 5.11, 5.11, 5.11, 5.11, 5.11, 5.11, 5.11, 5.12,5.12, 5.12, 5.12, 5.12, 5.12, 5.12, 5.12, 5.12, 5.12, 5.12, 5.12, 5.12, 5.12, 5.12,5.12, 5.12, 5.12, 5.13, 5.13, 5.13, 5.13, 5.13, 5.13, 5.13, 5.13, 5.13, 5.13, 5.14,5.14, 5.14, 5.14, 5.14, 5.14, 5.14, 5.14, 5.14, 5.14, 5.14, 5.14, 5.14, 5.15, 5.15,5.15, 5.15, 5.15, 5.15, 5.15, 5.16, 5.16, 5.16, 5.16, 5.16, 5.16, 5.17, 5.17, 5.17,5.17, 5.17, 5.17, 5.17, 5.17, 5.17, 5.18, 5.18, 5.18, 5.18, 5.18, 5.19, 5.19, 5.19,5.19, 5.19, 5.19, 5.19, 5.19, 5.20, 5.20, 5.20, 5.21, 5.21, 5.21, 5.22, 5.22, 5.22,5.22, 5.22, 5.22, 5.22, 5.22, 5.23, 5.23, 5.23, 5.23, 5.23, 5.24, 5.24, 5.24, 5.24,5.24, 5.24, 5.25, 5.25, 5.25, 5.26, 5.26, 5.27, 5.27, 5.27, 5.27, 5.27, 5.28, 5.28,5.29, 5.29, 5.30, 5.30, 5.31, 5.31, 5.32, 5.32, 5.32, 5.32, 5.33, 5.34, 5.34, 5.34,5.35, 5.35, 5.37, 5.38, 5.38, 5.40, 5.41, 5.46, 5.49, 5.50, 5.51, 5.51, 5.62, 5.64]:
XX := sort(X):

>    nops(X);
Percentile(X, 0.025), Percentile(X, 0.975);

1000

4.58, 5.30

>    xtics := [seq(4.25 + 0.1*k, k = 0 .. 14)]:
ytics := [seq(0.2*k, k = 1 .. 12)]:
B1 := HistogramFill(X,4.20 .. 5.70,15):
B2 := plot([[4.56,-0.06],[4.58,0], [4.6, -0.06]],x = 4.12 .. 5.77, y = -0.07 .. 2.56, xtickmarks=xtics, ytickmarks=ytics,labels=[``,``], color=black, thickness=2):
txt1 := textplot([4.58,-0.07,`4.58`],align=BELOW, font=[TIMES,BOLD,12]):
B3 := plot([[5.28,-0.06],[5.30,0], [5.32, -0.06]],x = 4.12 .. 5.77, y = -0.07 .. 2.56, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``], color=black, thickness=2):
txt2 := textplot([5.30,-0.07,`5.30`],align=BELOW, font=[TIMES,BOLD,12]):
display({B1, B2, B3,txt1,txt2});

[Maple Plot]

Figure 5.4-1(a)  Histogram of the Trimmed Means

>    F := ClassFreq(X,4.20 .. 5.70,15);

F := [1, 3, 10, 17, 55, 107, 200, 251, 204, 90, 38, 16, 4, 2, 2]

>    xbar := Mean(X);
var := Variance(X);

xbar := 4.943410000

var := .3087194194e-1

>    P := [seq(NormalP(0,1,k/1001), k = 1 .. 1000)]:

>    xtics := [seq(4.2 + 0.1*k, k = 0 .. 15)]:
#P1 := ScatPlotCirc(X,P):
P1 := ScatPlotPoint(X,P):
P2 := plot([[4.58,-0.1], [4.58, 0.1]], color=black,  xtickmarks=xtics, labels=[``,``], thickness=2):
txt1 := textplot([4.58,0.3,`4.58`],align=ABOVE, font=[TIMES,ROMAN,12]):
P3 := plot([[5.30,-0.1], [5.30, 0.1]], color=black, xtickmarks=xtics, labels=[``,``], thickness=2):
txt2 := textplot([5.30,0.3,`5.30`],align=ABOVE, font=[TIMES,ROMAN,12]):
display({P1, P2, P3,txt1,txt2});

[Maple Plot]

Figure 5.4-1(b) Quantiles of N (0,1) Versus Quantiles of Trimmed Means

>   

At this point you may return to the Menu or you may do a simulation similar to the one that was used to generate the figures above.

Return to Menu

Here are the original 40 observations of W .

>    W := [-7.34, -5.92, -2.98, 0.19, 0.77, 0.95, 2.86, 3.17, 3.76, 4.20, 4.20, 4.27, 4.31, 4.42, 4.60, 4.73, 4.84, 4.87, 4.90, 4.96, 4.98, 5.00, 5.09, 5.09, 5.14, 5.22, 5.23, 5.42, 5.50, 5.83, 5.94, 5.95, 6.00, 6.01, 6.24, 6.82, 9.62, 10.03, 18.27, 93.62]:

To find the value of the trimmed mean, we use the mean of the middle 10 order statistic which can be calculated as follows.

>    Y := [seq(W[m], m = 16 .. 25)];
Mean(Y);

Y := [4.73, 4.84, 4.87, 4.90, 4.96, 4.98, 5.00, 5.09, 5.09, 5.14]

4.960000000

To perform the simulation, create an empirical distribution, with probability of 1/40 on each observation.

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

empPDF := [-7.34, 1/40, -5.92, 1/40, -2.98, 1/40, .19, 1/40, .77, 1/40, .95, 1/40, 2.86, 1/40, 3.17, 1/40, 3.76, 1/40, 4.20, 1/40, 4.20, 1/40, 4.27, 1/40, 4.31, 1/40, 4.42, 1/40, 4.60, 1/40, 4.73, 1/40...
empPDF := [-7.34, 1/40, -5.92, 1/40, -2.98, 1/40, .19, 1/40, .77, 1/40, .95, 1/40, 2.86, 1/40, 3.17, 1/40, 3.76, 1/40, 4.20, 1/40, 4.20, 1/40, 4.27, 1/40, 4.31, 1/40, 4.42, 1/40, 4.60, 1/40, 4.73, 1/40...
empPDF := [-7.34, 1/40, -5.92, 1/40, -2.98, 1/40, .19, 1/40, .77, 1/40, .95, 1/40, 2.86, 1/40, 3.17, 1/40, 3.76, 1/40, 4.20, 1/40, 4.20, 1/40, 4.27, 1/40, 4.31, 1/40, 4.42, 1/40, 4.60, 1/40, 4.73, 1/40...
empPDF := [-7.34, 1/40, -5.92, 1/40, -2.98, 1/40, .19, 1/40, .77, 1/40, .95, 1/40, 2.86, 1/40, 3.17, 1/40, 3.76, 1/40, 4.20, 1/40, 4.20, 1/40, 4.27, 1/40, 4.31, 1/40, 4.42, 1/40, 4.60, 1/40, 4.73, 1/40...

Now generate 1000 resamples of size 40 and for each, calculate the value of the trimmed mean.

>    randomize:

>    for k from 1 to 1000 do
  XX := sort(DiscreteS(empPDF, 40)):
  Y := [seq(XX[m], m = 16 .. 25)]:
  XXbartrim[k] := Mean(Y):
od:

>    trimxbars := [seq(XXbartrim[k], k = 1 .. 1000)]:

>    Mean(trimxbars);

4.942494000

To find an approximate 95% confidence interval for theta   we use the 0.025 and 0.975 percentiles of the sample of trimmed x -bars (trimxbars).

>    lowerbound := Percentile(trimxbars, 0.025);
upperbound := Percentile(trimxbars, 0.975);

lowerbound := 4.537025000

upperbound := 5.303775000

>    lb := evalf(lowerbound,3);
ub := evalf(upperbound,3);

lb := 4.54

ub := 5.30

>    xtics := [seq(4.25 + 0.1*k, k = 0 .. 14)]:
ytics := [seq(0.2*k, k = 1 .. 14)]:
B1 := HistogramFill(trimxbars,4 .. 5.90,19):
B2 := plot([[lowerbound-.02,-0.06],[lowerbound,0], [lowerbound+.02, -0.06]],x = 4.12 .. 5.77, y = -0.07 .. 2.76, xtickmarks=xtics, ytickmarks=ytics,labels=[``,``], color=black, thickness=2):
txt1 := textplot([lb,-0.07,lb],align=BELOW, font=[TIMES,ROMAN,14]):
B3 := plot([[upperbound-.02,-0.06],[upperbound,0], [upperbound+.02, -0.06]],x = 4.12 .. 5.77, y = -0.07 .. 2.66, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``], color=black, thickness=2):
txt2 := textplot([ub,-0.07,ub],align=BELOW, font=[TIMES,ROMAN, 14]):
display({B1, B2, B3,txt1,txt2});

[Maple Plot]

LIKE: Figure 5.4-1(a)  Histogram of the Trimmed Means

>   

>    trimxbars := sort(trimxbars):

>    xtics := [seq(3.8 + 0.1*k, k = 0 .. 20)]:
P1 := ScatPlotPoint(trimxbars,P):
P2 := plot([[lowerbound,-0.1], [lowerbound, 0.1]], color=black,  xtickmarks=xtics, labels=[``,``], thickness=2):
txt1 := textplot([lb,0.3,lb],align=ABOVE, font=[TIMES,ROMAN,12]):
P3 := plot([[upperbound,-0.1], [upperbound, 0.1]], color=black, xtickmarks=xtics, labels=[``,``], thickness=2):
txt2 := textplot([ub,0.3,ub],align=ABOVE, font=[TIMES,ROMAN,12]):
display({P1, P2, P3,txt1,txt2});

[Maple Plot]

LIKE: Figure 5.4-1(b) Quantiles of N (0,1) Versus Quantiles of Trimmed Mean

>   

Return to Menu