proc iml;

start MC_PI(nsim, seed);
  temp_mat = J(nsim,2,0);
  call randseed(seed);
  call randgen(temp_mat,'uniform');
  temp_mat=temp_mat||(temp_mat[,2]<=sqrt(J(nsim,1,1)-
                        temp_mat[,1]##2));
 pi_over4 = temp_mat[+,3]/nsim;
pi_est = 4*pi_over4;
se_est = 4*sqrt(pi_over4*(1-pi_over4)/nsim);
pi_LCL = pi_est - 2*se_est;
pi_UCL = pi_est + 2*se_est;
*-----------------------------------------------------------;
print 'Estimating PI using MC integration method with' nsim
'data points';
print pi_est se_est pi_LCL pi_UCL;
finish MC_PI;


run MC_PI(400, 12345);
run MC_PI(1600, 12345);
run MC_PI(4000, 12345);

/*******/
quit;


/* example code with T/F indicator coerced to 1/0 */
data test;
  input x y;
  indicator = (y <= sqrt(1-x*x));
  if (y <= sqrt(1-x*x)) then statement=" TRUE";
  else statement="FALSE";
  datalines;
.5 .5
.75 1
;
run;
proc print data=test;
run;



/* alternate version of MC_PI module 
   - steps expanded 
*/
proc iml;

start MC_PI2(nsim, seed);
  temp_mat = J(nsim,2,0);
  call randseed(seed);
  call randgen(temp_mat,'uniform');

  ones_vector = J(nsim,1,1);
  x_sq = temp_mat[,1]##2;
  y = temp_mat[,2];
  under_circle = (y <=sqrt(ones_vector - x_sq));

  temp_mat=temp_mat|| under_circle;

 pi_over4 = temp_mat[+,3]/nsim;
 pi_est = 4*pi_over4;
 se_est = 4*sqrt(pi_over4*(1-pi_over4)/nsim);
 pi_LCL = pi_est - 2*se_est;
 pi_UCL = pi_est + 2*se_est;
*-----------------------------------------------------------;
 print 'Estimating PI using MC integration method with' \
       nsim 'data points';
 print pi_est se_est pi_LCL pi_UCL;
finish MC_PI2;