/* SAS program to accompany the article
   "Generate points uniformly inside a d-dimensional ball"
   Rick Wicklin, The DO Loop blog, published 06APR2016
   http://blogs.sas.com/content/iml/2016/04/06/generate-points-uniformly-in-ball.html

   This program simulates random uniform points in d-dimensional ball of radius R 
*/

/* DATA step to simulate random uniform points in d-dimensional ball of radius R */
%let dim = 3;
data Ball;
array X[&dim];
call streaminit(12345);
N = 1000;                            /* sample size */
radius = 2;                          /* use circle of radius 2 */
do i = 1 to N;
   do j = 1 to dim(X);
      X[j] = rand("Normal");       /* Y ~ MVN(0, I) */
   end;
   norm = sqrt(uss(of X[*]));
   u = rand("Uniform");           /* U ~ U(0,1) */
   r = radius * u**(1/&dim);               /* r proportional to d_th root of U */
   do j = 1 to dim(X);
      X[j] = r * X[j] / norm;
   end;
   output;
end;
run;

proc sgscatter data=Ball;
matrix X:;
run;

/* PROC IML program to simulate random uniform points in d-dimensional ball of radius R */
proc iml;
call randseed(12345);
d = 3;                         /* dimension = numbe of variables */
N = 1000;                      /* sample size = number of obs */
maxR = 2;                      /* radius of ball   */
Y = randfun(N // d, "Normal"); /* Y ~ MVN(0, I(d)) */
u = randfun(N, "Uniform");     /* U ~ U(0,1)       */
r = maxR * u##(1/d);           /* need density proportional to volume */
X = r # Y / sqrt(Y[,##]);      /* Y / ||Y|| is on d-sphere, scale by r */

/* if you want to construct an animated gif of the 3-D rotating cloud,
   use SAS/IML Studio and uncomment the following statements */

/*
* make sure bounds in each coordinate are from [-R, R];
fake = {-2 -2 .,
        -2 . -2,
         2  2 .,
         2  . 2};
X = X // fake;

declare RotatingPlot plot = 
     RotatingPlot.Create("unifB3", X[,1], X[,2], X[,3]);
plot.SetTitleText("Random Uniform Points in 3-Ball", true);
plot.SetGraphAreaBackgroundColor(WHITE);
plot.SetWindowPosition(50, 0, 50, 100);
plot.SetAxisTickRange( XAXIS, -2, 2 );
plot.SetAxisTickRange( YAXIS, -2, 2 );
plot.SetAxisTickRange( ZAXIS, -2, 2 );
plot.SetAxesLocation( MINIMA );
pause "CTRL+B to add bounding box, then RESUME";

* save images to disk for animation;
path = 'C:\Temp\AnimGif\Sim3D\';               * directory to store images;
plot.SaveToFile(path+"Sim000.bmp");            * save the initial configuration;
do i = 1 to 63;    
   plot.Rotate( 3.14159/32, 0, 1, 1 ); 
   suffix = trim(putn(i, "Z3."));              * 001, 002, ..., 063;
   plot.SaveToFile(path+"Sim"+suffix+".bmp");  * save image to Sim001.bmp;
end;

* You can then use a free public utility like GiftedMotion to create the animated gif;
*/