SIMSYSTEM Procedure

(View the complete code for this example.)

This example demonstrates how to use the table of parameters that PROC SIMSYSTEM displays to reconstruct in DATA step code random variates that have the specified moments for each distribution. Starting with the Pearson system, the moment combinations that are specified for PROC SIMSYSTEM in the following code are chosen to span the families of distributions in the Pearson family, plus the normal distribution, as the parameters table (Output 31.5.2) that PROC SIMSYSTEM produces for these moment combinations confirms:

%let Mean   = -4.5 -3.5 -2.5 -1.5 -0.5  0.5        1.5  2.5        3.5  4.5;
%let Stddev =  1    2    3    4    5    6          7    8          9   10;
%let Skew   =  0   -1.5  1.8 -1.3  0   -2.00831604 0.1 -1.04235831 1.7  0;
%let Kurt   =  3    5.7  7.4  3    2.1  9.05       6    5.15       8.6  3.2;
proc simsystem system=pearson;
   moments mean=&Mean StdDev=&Stddev skew=&Skew kurt=&Kurt;
   ods output Parameters=Parm;
run;

The ODS OUTPUT statement in the preceding PROC SIMSYSTEM code saves the parameters for this sequence of distributions in a data set named Parm. The following PROC PRINT code produces the selected contents of this data set shown in Output 31.5.3. As you can see, the exact numeric values of the parameters that describe each distribution in the Notes column are available in the Shap1Val and Shap2Val variables.

proc print data=Parm noobs;
   var Family Notes Shap: Shift Scale;
run;

As discussed in Parameter Values for Pearson Distributions, you can use the values of the Shap1Val and Shap2Val variables in DATA step code with the formulas given in Table 3 to compute random samples that have the specified moments, as shown in the following DATA step code:

%let N = 1000000;
data mylib.Parm; set Parm; run;
proc ds2 sessref=mysess;
data Sim / overwrite=yes;
   dcl int i;
   dcl double y;
   method run();
      set Parm;
      streaminit(15531);
      do i=1 to &N;
         select (trim(left(Family)));
            when ('Type I'   )
            when ('Type I(J)')
            when ('Type I(U)') y =   rand('beta'  ,Shap1Val,Shap2Val);
            when ('Type II'  ) y =   rand('beta'  ,Shap1Val,Shap1Val);
            when ('Type III' ) y =   rand('gamma' ,Shap1Val);
            when ('Type IV'  ) y = .;
            when ('Type V'   ) y = 1/rand('gamma' ,Shap1Val);
            when ('Type VI'  ) y =   rand('f'     ,Shap1Val,Shap2Val);
            when ('Type VII' ) y =   rand('t'     ,Shap1Val);
            when ('Normal'   ) y =   rand('normal');
         end;
         y = Scale*y + Shift;
         output;
      end;
   end;
enddata;
run;
quit;

The following code computes and displays (in Output 31.5.4) the first four moments for each distribution sampled by the preceding PROC DS2 code. These moments match the specified ones reasonably well.

proc mdsummary data=mylib.Sim;
   var y;
   groupby Family Notes Mean StdDev Skewness Kurtosis;
   output out=mylib.Sim_summary;
run;

proc cas;
   session mysess;
   table.alterTable /
      name="Sim_summary"
      columns={
         {name="_Mean_"       rename="sMean"},
         {name="_Std_"        rename="sStdDev"},
         {name="_Skewness_"   rename="sSkewness"},
         {name="_Kurtosis_"   rename="sKurtosis"}
      };
   table.update /
      table="Sim_summary"
      set={{var="sKurtosis" value="sKurtosis + 3"}};
quit;

proc sort data=mylib.Sim_summary out=Sim_summary; by Family; run;
proc print data=Sim_summary noobs;
   var Family Mean sMean StdDev sStdDev Skewness sSkewness Kurtosis sKurtosis;
   format _NUMERIC_ 6.3;
run;

Finally, for the Johnson system, the following code performs steps precisely like those shown earlier for the Pearson system, except that the moment combinations are chosen to represent the families of the Johnson system, as shown in Output 31.5.5 and Output 31.5.6:

%let Mean   =  -2   -1    0            1  2  ;
%let Stddev =   1    2    3            4  5  ;
%let Skew   =  -1.6  2.5 -1.2922847983 0 -0.7;
%let Kurt   =   7.5  8.4  6.1093660061 3  7  ;
proc simsystem system=johnson;
   moments mean=&Mean StdDev=&Stddev skew=&Skew kurt=&Kurt;
   ods output Parameters=Parm;
run;

Turning to the Johnson system, as discussed in Parameter Values for Johnson Distributions, you can use the values of the Gamma and Delta columns of the "Parameters for Johnson Distributions" table in DATA step code with the formulas given in Table 5 to compute random samples that have the specified moments, as shown in the following code:

%let N = 1000000;
data mylib.Parm; set Parm; run;
proc ds2 sessref=mysess;
data Sim / overwrite=yes;
   dcl int i;
   dcl double y z;
   method run();
      set Parm;
      streaminit(15531);
      do i = 1 to &N;
         z = rand('normal');
         select (trim(left(Family)));
            when ('SB(1)')
            when ('SB(2)') y = 1/(1 + exp(-(z - Gamma)/Delta));
            when ('SL'   ) y =        exp(  z         /Delta) ;
            when ('SU'   ) y =       sinh( (z - Gamma)/Delta) ;
            when ('SN'   ) y =              z                 ;
         end;
         y = Scale*y + Shift;
         output;
      end;
   end;
enddata;
run;
quit;

proc mdsummary data=mylib.Sim;
   var y;
   groupby Family Mean StdDev Skewness Kurtosis;
   output out=mylib.Sim_summary;
run;

proc cas;
   session mysess;
   table.alterTable /
      name="Sim_summary"
      columns={
         {name="_Mean_"       rename="sMean"},
         {name="_Std_"        rename="sStdDev"},
         {name="_Skewness_"   rename="sSkewness"},
         {name="_Kurtosis_"   rename="sKurtosis"}
      };
   table.update /
      table="Sim_summary"
      set={{var="sKurtosis" value="sKurtosis + 3"}};
quit;

proc sort data=mylib.Sim_summary out=Sim_summary; by Family; run;
proc print data=Sim_summary noobs;
   var Family Mean sMean StdDev sStdDev Skewness sSkewness Kurtosis sKurtosis;
   format _NUMERIC_ 6.3;
run;

As was the case with the Pearson system shown earlier, the moments for the samples that the PROC DS2 code computes, shown in Output 31.5.7, match the specified moments reasonably well.