*! version 1.0.0 14may1995 program define sdist version 4.0 preserve /* Save all the current stata stuff */ clear /* Clear everything so no conflicts */ quietly { set obs 501 gen x=. gen y=. gen xbar1=. gen xbar2=. gen s21=. gen s22=. gen stat=.} global D_sm11 " This lab shows that we can use the percentiles of the Z, t, Chi-square, and F curves to determine a range of likely values for a transformed statistic. Later we will use these likely values to do statistical inferences (confidence intervals and tests of hypotheses) about parameters." wdctl static D_sm11 5 5 295 32 global D_sm12 " The lab does what is called a Monte Carlo study, that is, it simulates 500 samples (or pairs of samples for the two sample t or F statistics) from a population (or pair of populations) having one of three possible distributions. For each sample or pair of samples, the lab calculates the value of a transformed statistic. It then draws a histogram of these 500 values and draws the theoretical distribution on top of the histogram. Finally the percentiles of the theoretical curve and the histogram are written on the graph so you can see how closely they agree." wdctl static D_sm12 5 40 295 56 global D_sm13 " You choose 1) Which statistic to study, 2) Which population to sample from, and 3) The sample size (or sizes if there are pairs of samples)." wdctl static D_sm13 5 100 295 16 global D_sm1 "n_1" global HJN_V1=5 wdctl static D_sm1 80 130 12 10 wdctl edit HJN_V1 95 130 16 10 global D_sm2 "n_2" global HJN_V2=5 wdctl static D_sm2 120 130 12 10 wdctl edit HJN_V2 135 130 16 10 global D_sm6 "Sample from" global D_sm7 "Uniform(0,1) Exponential(1) Normal(0,1)" wdctl static D_sm6 160 120 60 10 global D_var1 "Normal(0,1)" wdctl ssimple D_var1 D_sm7 160 130 55 48 global D_sm8 "Statistic" global D_sm9 "Z One-sample-t Two-sample-t Chi-square F" wdctl static D_sm8 230 120 45 10 global D_var2 "Z" wdctl ssimple D_var2 D_sm9 230 130 55 60 wdctl button "Run" 5 160 30 14 D_b1 wdctl button "Close" 40 160 30 14 D_b2 wdctl button "Help" 75 160 30 14 D_b3 help global D_b1 "sdistdr" global D_b2 "exit 1234" global D_b3 "whelp sdist" cap noi wdlg "Sampling Distributions" $D_dlgx $D_dlgy 320 240 restore end program define sdistdr version 4.0 local n1=$HJN_V1 local n2=$HJN_V2 scalar nsamps=1 if "$D_var2"=="F" | "$D_var2"=="Two-sample-t" {scalar nsamps=2} gph open gph pen 2 if "$D_var1"=="Uniform(0,1)" { scalar sig2 = 1.0 / 12. scalar mu = 0.5} if "$D_var1"=="Exponential(1)" { scalar sig2 = 1.0 scalar mu = 1.0} if "$D_var1"=="Normal(0,1)" { scalar sig2 = 1.0 scalar mu = 0.0} scalar sig = sqrt(sig2) /* Start while loop */ replace xbar1=0 replace s21=0 if nsamps==1 { gph text 1000 16000 0 0 Generating 500 Samples From a $D_var1 Population gph text 2000 16000 0 0 and Calculating $D_var2 for Each Sample} else { gph text 1000 16000 0 0 Generating 500 Pairs of Samples From a $D_var1 Population gph text 2000 16000 0 0 and Calculating $D_var2 for Each Pair of Samples} gph box 6000 10000 8000 20000 4 scalar i = 1 while i <= `n1'{ local mi=i local bx2=int(10000 + 10000*(i/`n1')) gph box 6000 10000 8000 `bx2' 1 if "$D_var1"=="Uniform(0,1)" { replace x = uniform()} if "$D_var1"=="Exponential(1)" { replace x = -log(uniform())} if "$D_var1"=="Normal(0,1)" { replace x = invnorm(uniform())} replace xbar1=xbar1+x replace s21=s21+x^2 scalar i=i+1 } replace s21=(s21-(xbar1^2/`n1'))/(`n1'-1) replace xbar1=xbar1/`n1' if nsamps==2 { replace xbar2=0 replace s22=0 gph box 6000 10000 8000 20000 4 scalar i = 1 while i <= `n2'{ local bx2=int(10000 + 10000*(i/`n2')) gph box 6000 10000 8000 `bx2' 1 if "$D_var1"=="Uniform(0,1)" { replace x = uniform()} if "$D_var1"=="Exponential(1)" { replace x = -log(uniform())} if "$D_var1"=="Normal(0,1)" { replace x = invnorm(uniform())} replace xbar2=xbar2+x replace s22=s22+x^2 scalar i=i+1 } replace s22=(s22-(xbar2^2/`n2'))/(`n2'-1) replace xbar2=xbar2/`n2' } if "$D_var2"=="Z" {replace stat = (xbar1-mu)/sqrt(sig2/`n1')} if "$D_var2"=="One-sample-t" {replace stat = (xbar1-mu)/sqrt(s21/`n1')} if "$D_var2"=="Chi-square" {replace stat = ((`n1'-1)*s21)/sig2} if "$D_var2"=="Two-sample-t" { scalar s2p=((`n1'-1)*s21+(`n2'-1)*s22)/(`n1'+`n2'-2) replace stat = (xbar1-xbar2)/sqrt(s2p*((1/`n1')+(1/`n2')))} if "$D_var2"=="F" {replace stat = s21/s22} set graphics off graph stat, bin(20) xlabel ylabel l1(" ") l2(" ") b1(" ") b2(" ") set graphics on local ymin = _result(1) local ymax = _result(2) local xmin = _result(3) local xmax = _result(4) replace x = `xmin' + ((_n-1)/500)*(`xmax'-`xmin') local ymax=`ymax'/((`xmax'-`xmin')/20) replace y = `ymax'*uniform() gph close gph open gph pen 1 graph y x, s(i) xlab ylab l1(" ") l2(" ") b1(" ") b2(" ") bbox(3000,0,23000,25000,850,400,0) gphconv if "$D_var2"=="Chi-square" { local df=`n1'-1 pdf_chi2 x y `df' local z1: display %5.2f invnchi(`df',0,.025) local z2: display %5.2f invnchi(`df',0,.050) local z3: display %5.2f invnchi(`df',0,.100) local z4: display %5.2f invnchi(`df',0,.250) local z5: display %5.2f invnchi(`df',0,.500) local z6: display %5.2f invnchi(`df',0,.750) local z7: display %5.2f invnchi(`df',0,.900) local z8: display %5.2f invnchi(`df',0,.950) local z9: display %5.2f invnchi(`df',0,.975)} if "$D_var2"=="Z" { replace y = exp(-(x^2)/2)/sqrt(2*_pi) local z1: display %5.2f invnorm(.025) local z9: display %5.2f -`z1' local z2: display %5.2f invnorm(.050) local z8: display %5.2f -`z2' local z3: display %5.2f invnorm(.100) local z7: display %5.2f -`z3' local z4: display %5.2f invnorm(.250) local z6: display %5.2f -`z4' local z5: display %5.2f 0} if "$D_var2"=="One-sample-t" | "$D_var2"=="Two-sample-t" { local df=`n1'-1 if"$D_var2"=="Two-sample-t" {local df=`n1'+`n2'-2} pdf_t x y `df' local z5: display %5.2f 0 local z9: display %5.2f invt(`df',.95) local z1: display %5.2f -`z9' local z8: display %5.2f invt(`df',.90) local z2: display %5.2f -`z8' local z7: display %5.2f invt(`df',.80) local z3: display %5.2f -`z7' local z6: display %5.2f invt(`df',.50) local z4: display %5.2f -`z6'} if"$D_var2"=="F" { local df1 = `n1'-1 local df2 = `n2'-1 pdf_f x y `df1' `df2' local z1: display %5.2f invfprob(`df1',`df2',.975) local z2: display %5.2f invfprob(`df1',`df2',.950) local z3: display %5.2f invfprob(`df1',`df2',.900) local z4: display %5.2f invfprob(`df1',`df2',.750) local z5: display %5.2f invfprob(`df1',`df2',.500) local z6: display %5.2f invfprob(`df1',`df2',.250) local z7: display %5.2f invfprob(`df1',`df2',.100) local z8: display %5.2f invfprob(`df1',`df2',.050) local z9: display %5.2f invfprob(`df1',`df2',.025)} * replace y = y*(`xmax'-`xmin')/20 * summarize y * local yy1 = _result(6) * if(`yy1' > `ymax') {local ymax=`yy1'} histgm stat `xmin' `xmax' 20 gph pen 3 lines x y gph text 1000 16000 0 0 Theoretical and Monte Carlo Distribution of $D_var2 if nsamps==1 {gph text 2000 16000 0 0 Sampling From $D_var1, n=`n1'} else {gph text 2000 16000 0 0 Sampling From $D_var1, n1=`n1', n2=`n2'} * gph text 21000 12000 0 -1 95% Fall in (`z1',`z2') : Curve sort stat local o1: display %5.2f (stat[13]+stat[14])/2 local o2: display %5.2f stat[26] local o3: display %5.2f stat[51] local o4: display %5.2f stat[126] local o5: display %5.2f stat[251] local o6: display %5.2f stat[376] local o7: display %5.2f stat[451] local o8: display %5.2f stat[476] local o9: display %5.2f (stat[488]+stat[489])/2 gph pen 2 * gph text 22000 12000 0 -1 95% Fall in (`z1',`z2') : Samples gph text 5000 26000 0 0 Percentiles of curve gph text 6000 26000 0 0 and histogram count if stat < `z1' local w1: display %4.3f _result(1)/501 count if stat < `z2' local w2: display %4.3f _result(1)/501 count if stat < `z3' local w3: display %4.3f _result(1)/501 count if stat < `z4' local w4: display %4.3f _result(1)/501 count if stat < `z5' local w5: display %4.3f _result(1)/501 count if stat < `z6' local w6: display %4.3f _result(1)/501 count if stat < `z7' local w7: display %4.3f _result(1)/501 count if stat < `z8' local w8: display %4.3f _result(1)/501 count if stat < `z9' local w9: display %4.3f _result(1)/501 gph text 8000 21500 0 -1 .025 `z1' `o1' `w1' gph text 9000 21500 0 -1 .050 `z2' `o2' `w2' gph text 10000 21500 0 -1 .100 `z3' `o3' `w3' gph text 11000 21500 0 -1 .250 `z4' `o4' `w4' gph text 12000 21500 0 -1 .500 `z5' `o5' `w5' gph text 13000 21500 0 -1 .750 `z6' `o6' `w6' gph text 14000 21500 0 -1 .900 `z7' `o7' `w7' gph text 15000 21500 0 -1 .950 `z8' `o8' `w8' gph text 16000 21500 0 -1 .975 `z9' `o9' `w9' gph line 4000 21000 4000 31000 1 gph line 4000 31000 17000 31000 1 gph line 17000 31000 17000 21000 1 gph line 17000 21000 4000 21000 1 gph close end exit