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

Warning, the name changecoords has been redefined

Elliot Tanis
April 16, 2001
T-distribution -- Equal Variances--  (Tev.mws)

First simulate N  samples of size n  from a N (0, VarX) distribution and N   samples of size m  from a N (0, VarY) distribution. Calculate the value of  " T" ,

t = (xbar-ybar-(0-0))/(sqrt(((n-1)*s[x]^2+(m-1)*s[y]^2)/(n+m-2))*sqrt(1/n+1/m))

and see if T  actually has a T -distribution with n  + m  - 2 degrees of freedom.  Also calculate the value of   " Z" ,

z = (xbar-ybar-(0-0))/sqrt(s[x]^2/n+s[y]^2/m)  .

How do the sizes of n  and m  affect the distribution of T ? How important is it for the distribution variances to be equal for T  to actually have a t-distribution with n  + m  - 2 degrees of freedom?

For this simulation the variances are equal.

>    randomize():

>    N := 500;
n := 6;
m := 18;
VarX := 16;
VarY := 16;

N := 500

n := 6

m := 18

VarX := 16

VarY := 16

>    for k from 1 to N do
X := NormalS(0, VarX, n):
Y := NormalS(0, VarY, m):
TT[k] := evalf((Mean(X) - Mean(Y))/sqrt(((n-1)*Variance(X) + (m-1)*Variance(Y))/(n + m - 2)*(1/n + 1/m))):
ZZ[k] := evalf((Mean(X) - Mean(Y))/sqrt(Variance(X)/n + Variance(Y)/m)):
od:
T := [seq(TT[j], j = 1 .. N)]:
Z := [seq(ZZ[j], j = 1 .. N)]:

>    MedianofTs := Median(T);
MedianofZs := Median(Z);
MeanofTs := Mean(T);
MeanofZs := Mean(Z);
SDofTs := StDev(T);
SDofZs := StDev(Z);
Min(T), Max(T);
Min(Z), Max(Z);

MedianofTs := .6055005461e-2

MedianofZs := .6191578299e-2

MeanofTs := .2103637932e-2

MeanofZs := .3075573076e-2

SDofTs := .9826343078

SDofZs := 1.078610371

-3.326584220, 3.613884909

-3.834562371, 3.173911686

>    display({BoxPlot(T, Z)}, title = `Box Plots for Z (above) and T (below)`, titlefont=[TIMES,BOLD,16]);

[Maple Plot]

>    xtics := [seq(-3.5 + k*0.5, k = 0 .. 14)]:
ytics := [seq(0.1*k, k = 1 .. 11)]:

>    P1 := HistogramFill(T, -3.5 .. 3.5, 28):
P2 := plot(TPDF(n+m-2, x), x = -3.5 .. 3.5, y = 0 .. 0.45, color=black, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]):
display({P1, P2}, title = `Histogram of T Observations\nwith T(n+m-2) Superimposed`, titlefont=[TIMES,BOLD,16]);

`WARNING: There are`, 1, ` data points out of plot range.`

[Maple Plot]

>    n+m-2;

22

>    Percentile(T, 0.05), Percentile(T, 0.95);
TP(n+m-2,0.05), TP(n+m-2,0.95);

-1.585194784, 1.572909827

-1.717144374, 1.717144374

>    P3 := HistogramFill(Z, -3.5 .. 3.5, 28):
P4 := plot(NormalPDF(0, 1, x), x = -3.5 .. 3.5, y = 0 .. 0.45, color=black, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]):
display({P3, P4}, title = `Histogram of Z Observations\nwith N(0, 1) p.d.f. Superimposed`, titlefont=[TIMES,BOLD,16]);

`WARNING: There are`, 2, ` data points out of plot range.`

[Maple Plot]

>    Percentile(Z, 0.05), Percentile(Z, 0.95);
NormalP(0, 1, 0.05), NormalP(0, 1, 0.95);

-1.730431809, 1.724857765

-1.644853627, 1.644853627

>    TQuants := [seq(TP(n+m-2, k/(N+1)), k = 1 .. N)]:

>    T := sort(T):
P5 := ScatPlotCirc(T, TQuants):
P6 := plot([[-3.5, -3.5], [3.5, 3.5]], x = -3.5 .. 3.5, y = -3.5 .. 3.5, color=black, thickness=2, xtickmarks=xtics, ytickmarks=xtics, labels=[``,``]):
display({P5, P6}, title=`q-q Plot of T(n+m-2) Quantiles\nVersus the T Order Statistics`, titlefont=[TIMES,BOLD,16]);

[Maple Plot]

>    n+m-2;

22

>    ZQuants := [seq(NormalP(0, 1, k/(N+1)), k = 1 .. N)]:

>    Z := sort(Z):
P7 := ScatPlotCirc(Z, ZQuants):
P6 := plot([[-3.5, -3.5], [3.5, 3.5]], x = -3.5 .. 3.5, y = -3.5 .. 3.5, color=black, thickness=2, xtickmarks=xtics, ytickmarks=xtics, labels=[``,``]):
display({P7, P6}, title=`q-q Plot of N(0, 1) Quantiles\nVersus the Z Order Statistics`, titlefont=[TIMES,BOLD,16]);

[Maple Plot]

>   

Welch showed that U = Z,  where

U = (Xbar-Ybar-(0-0))/sqrt(S[x]^2/n+S[y]^2/m)  ,

really has an approximate T  distribution with the number of degrees of freedom given by the greatest integer (or the floor) in the following expression:

v = (s[x]^2/n+s[y]^2/m)^2/((s[x]^2/n)^2/(n-1)+(s[y]^2/m)^2/(m-1))

We shall now fit a t -distribution to the Z  data using as the number of degrees of freedom Welch's suggestion but with distribution variances rather than sample variances.

>    v := floor((VarX/n + VarY/m)^2/((VarX/n)^2/(n-1) + (VarY/m)^2/(m-1)));

v := 8

>    P1 := HistogramFill(Z, -3.5 .. 3.5, 28):
P2 := plot(TPDF(v, x), x = -3.5 .. 3.5, y = 0 .. 0.45, color=black, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]):
display({P1, P2}, title=`Histogram of U Observations\nwith T(v) Superimposed`, titlefont=[TIMES,BOLD,16]);

`WARNING: There are`, 2, ` data points out of plot range.`

[Maple Plot]

>    v;

8

>    Percentile(Z, 0.05), Percentile(Z, 0.95);
TP(v, 0.05), TP(v, 0.95);

-1.730431809, 1.724857765

-1.859548038, 1.859548038

>    TvQuants := [seq(TP(v, k/(N+1)), k = 1 .. N)]:

>    Z := sort(Z):
P5 := ScatPlotCirc(Z, TvQuants):
P6 := plot([[-3.5, -3.5], [3.5, 3.5]], x = -3.5 .. 3.5, y = -3.5 .. 3.5, color=black, thickness=2, xtickmarks=xtics, ytickmarks=xtics, labels=[``,``]):
display({P5, P6}, title=`q-q Plot of T(v) Quantiles\nVersus the U Order Statistics`, titlefont=[TIMES,BOLD,16]);

[Maple Plot]

>    v;

8

Return to Menu

>