/*****************************************************************/
/*                                                               */
/*  random-means.jsl is a script that simulates the following    */
/*  process:                                                     */
/*                                                               */
/*     Draw a sample of size n using SRS from a population       */
/*       data set with N rows called Y (stored in the current    */
/*       data table, and compute the mean of these n draws;      */
/*     repeat this process M times and store these M means       */ 
/*       in a column called Mean in a data table called          */
/*       Results.                                                */
/*                                                               */
/*  In the language of the populations-and-samples diagrams      */
/*  from class, Y is the population and Mean is the imaginary    */
/*  data set.                                                    */
/*                                                               */
/*  When the population size N is a lot larger than the sample   */
/*  size n, SRS is approximately the same as IID sampling, so    */
/*  we can use the formulas in class for the expected value and  */
/*  standard error of a mean under IID sampling to get an idea   */
/*  of what the mean and SD of the column called Mean should be, */
/*  and if n is large the distribution of values in Mean should  */
/*  be close to normal by the Central Limit Theorem.             */
/*                                                               */
/*****************************************************************/

// Reference data table;

dt=current data table();

dt<<ClearRowStates;

// Specify sample size

n=5; // size of sample drawn from population;

// Specify the number of samples;

M=500; // number of simulated means in imaginary data set;

// Preparing a new data table for the results to be stored;

ResultsDT=NewTable("Results.jmp",
                     NewColumn("Mean"),
                   );

For(i=1,i<=M,i++,
     currentdatatable(dt);
     eval(substitute(expr(current data table()<<SelectRandomly(ss)), 
       expr(ss), n));

     dt<<Invert Row Selection;
     dt<<Hide;
     dt<<Exclude;
     Dist=Distribution(Continuous Distribution(Column(:Y)));

     RDist=Dist<<report;

     stats=RDist["Moments"][NumberColBox(1)]<<get as matrix;
       
     ResultsDT<<AddRow(1);
       Column(ResultsDT,1)[row()] = stats[1,1]; // Mean;
     Dist<<Close Window;
     dt<<ClearRowStates;
   );