Last week I showed how to generate random points uniformly inside a 2-d circular region.
That article showed that the distance of a point to the circle's center cannot be distributed uniformly. Instead, you should use the square root of a uniform variate to generate 2-D distances to the origin.
This result generalizes to higher dimensions. If you want to simulate points uniformly in the *d*-dimensional ball, generate the radius to be proportional to the *d*th-root of a uniform random variate.

Although it is possible to generalize the polar mapping and obtain a parameterization of the *d*-dimensional ball, there is an easier approach to simulate points in a ball. The basic idea is to simulate *d* independent standardized normal variables, project them radially onto the unit sphere, and then adjust their distance to the origin appropriately. You can find this algorithm in many textbooks (Devroye, 1986; Fishman, 1996), but
Harman and Lacko (2010) summarize the process nicely. To paraphrase:

If Y is drawn from the uncorrelated multivariate normal distribution, then S = Y / ||Y|| has the uniform distribution on the unitd-sphere. Multiplying S by U^{1/d}, where U has the uniform distribution on the unit interval (0,1), creates the uniform distribution in the unitd-dimensional ball.

It is easy to carry out this operation in SAS. In Base SAS, you can write a DATA step that uses arrays. The following program provides a SAS/IML solution. I used the programming features of SAS/IML Studio to create a 3-D rotating scatter plot of the simulated cloud of points. You can download the SAS program that contains the DATA step and the IMLPlus program that creates the animation.

proc iml; call randseed(12345); d = 3; /* dimension = number of variables */ N = 1000; /* sample size = number of obs */ radius = 2; /* radius of circle */ Y = randfun(N // d, "Normal"); /* Y ~ MVN(0, I(d)) */ u = randfun(N, "Uniform"); /* U ~ U(0,1) */ r = radius * u##(1/d); /* r proportional to d_th root of U */ X = r # Y / sqrt(Y[,##]); /* Y[,##] is sum of squares for each row */ /* X contains N random uniformly distributed points in the d-ball */ |

## 2 Comments

Pingback: Rotation matrices and 3-D data - The DO Loop

Great algorithm, Simple and elegant