/* program to perform some sampling experiments                                  c. flinn,  10/99 */

output file = d:\r\class\samp_1.our reset;

/* let the population distribution be normal with mean 5 and variance 25 */

/* for sample size n = 2, draw 1000 samples */

x = 5 + 5*rndn(5,1000);

@ compute the mean of each of the 1000 samples @
xbar = meanc(x);
@plot the histogram of xbar @
library pgraph;
graphset;
_pdate = "";
xtics(-5,15,5,0);
_pbarwid = 1;
title("Distribution of Sample Means\LNormal Population\LN = 5");
histp(xbar,100);
@ repeat for a sample size of 20 @
x = 5 + 5*rndn(20,1000);
xbar = meanc(x);
title("Distribution of Sample Means\LNormal Population\LN = 20");
histp(xbar,100);

/* repeat experiment with exponentially distributed random variable with same
   mean and variance as normal r.v. above.
   If lambda = .2, then mean(x) is 5 and var(x) is 25  */

/* to generate exponential r.v.s, use uniform random number generator 
 
    p = 1 - exp(-lambda*x)   is exponential c.d.f. 
    draw p from uniform on [0,1], invert to get value of x, or
    x = - ln(1-p)/lambda                    */

lambda = .2;
fn u_to_e(u_rv) = - ln(1-u_rv) / lambda;

@ experiment with N = 5 @
x = u_to_e(rndu(5,1000));
xbar = meanc(x);
title("Distribution of Sample Means\LExponential Population\LN=5");
histp(xbar,100);
@ repeat for N = 20 @
x = u_to_e(rndu(20,1000));
xbar = meanc(x);
title("Distribution of Sample Means\LExponential Population\LN=20");
histp(xbar,100);

/* we can use the computer to numerically evaluate sampling distributions 
    consider the exponential population example  with N = 5*/
@ draw a very large numbers of samples of size 5, and then enumerate the proportion
     with a mean value less than or equal to a large number of points of evaluation @
xx = u_to_e(rndu(5,100000));  /* 	WARNING: DON'T TRY THIS AT HOME 
                                                         OR EVEN IN OUR COMPUTER LAB!!! */
xxbar = meanc(xx);
p_eval = zeros(2001,2);
p_eval[.,1] = seqa(0,.01,2001);
i = 2;
do until i gt 2001;
      p_eval[i,2] = meanc( p_eval[i-1,1] .<= xxbar .and p_eval[i,1] .> xxbar);
      i = i+1;
endo;

title("Estimated p.d.f. of sampling distribution of mean\LExponential Pop\lLambda = .2");
xy(p_eval[.,1],p_eval[.,2]); 


end;