/* 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; */