/*****************************************************************************************************************

SAS file name: Central_Limit_Theorem.sas
File location: 
_________________________________________________________________________________________________________________

Purpose: To exaplin and visualize the Central Limit Theorem using SAS data step code, IML and PROC SGPLOT
Author: Peter Clemmensen
Creation Date: 01/08/2017

This program supports the blog post "Visualize the Central Limit Theorem in SAS" on SASnrd.com

*****************************************************************************************************************/

/* Tabulated distribution example */

/* Draw samples from unfair die rolls */
%let NumSample=100;
%let SampleSize=10000;
data DieRolls; 
   call streaminit(321);
   do sample=1 to &NumSample;
      do n=1 to &SampleSize;
          x=rand("Table", 0.2, 0.1, 0.1, 0.3, 0.2, 0.1);
          output;
      end;
   end;
run;

/* Visualize distribution of die rolls */
title 'Die Roll Frequency Distribution';
proc sgplot data=DieRolls;
   vbar x;
   xaxis values=(1 to 6);
run;
title;

/* Calculate sample means */
proc means data=DieRolls noprint;
   class sample;
   output out=DieRollMeans(where=(_STAT_='MEAN' and _TYPE_=1));
run;

/* Visualize samlping distribution of the mean */
title 'Sampling Distribution of the Mean';
proc sgplot data=DieRollMeans noautolegend;
   histogram x / scale=count;
   density x / type=normal;
run;
title;

/* Poisson distribution example */

/* Draw samples from Poisson distribution with parameter lambda */
%let NumSample=100;
%let SampleSize=10000;
%let lambda=3;
data PoissonSamples; 
   call streaminit(321);
   do sample=1 to &NumSample;
      do n=1 to &SampleSize;
          po=rand("Poisson", &lambda);
          output;
      end;
   end;
run;

/* Visualize simulaed density of poisson values */
title "Simulated Poisson distribution with (*ESC*){unicode lambda}=&lambda";
proc sgplot data=PoissonSamples;
   vbar po / barwidth=1;
   xaxis display=(nolabel);
run;
title;

/* Calculate sample means */
proc means data=PoissonSamples noprint;
   class sample;
   output out=PoissonSampleMeans(where=(_STAT_='MEAN' and _TYPE_=1));
run;

/* Visualize samlping distribution of the mean */
title 'Sampling distribution of the mean';
proc sgplot data=PoissonSampleMeans noautolegend;
   histogram po / scale=count;
   density po / type=normal;
   xaxis display=(nolabel);
run;
title;