15
$\begingroup$

A scheme to generate random variates distributed uniformly in $S^2=\{x\in \mathbb{R}^n \mid \|x\|_2=1\}$ is well known: generate a standard normal variate in $\mathbb{R}^n$ and normalize it to unit norm.

Is there a similarly simple and clever procedure to simulate uniformly distributed variates on the $1$-ball $S^1=\{x \in \mathbb{R}^n \mid \|x\|_1=1\}$?

|cite|improve this question
$\endgroup$
  • 1
    $\begingroup$ Have you tried the same procedure but normalizing using the 1-norm? $\endgroup$ – Raskolnikov Mar 6 '11 at 14:59
  • 1
    $\begingroup$ I fear that that might not work. The normal-variate approach generates points uniformly according to their 'angle'; in $R^2$, for example, their angles are distributed uniformly for $\theta \in [0,2\pi[.$ However, a small line element of the 1-ball in $R^2$ doesn't work that way; near a 'corner', it has $\mathrm{d}l = \frac{\sqrt{2}}{2} \mathrm{d}\theta,$ whereas in the middle of a line segment, we have $\mathrm{d}l = \mathrm{d}\theta.$ $\endgroup$ – Gerben Mar 6 '11 at 15:20
  • 3
    $\begingroup$ To make sure we're talking about the same thing here (since you didn't specify any measure): by "uniformly distributed", do you mean that the probability of obtaining $x$ within some subset of $S^1$ should be proportional to the surface (hyper)area of that subset as measured using the standard metric in $\mathbb{R}^n$? $\endgroup$ – joriki Mar 6 '11 at 15:41
  • $\begingroup$ [I'm sorry, but the factor of 1/2 (below the \sqrt) in my response should not be there, yet I can't edit my post anymore.] $\endgroup$ – Gerben Mar 6 '11 at 15:50
  • 4
    $\begingroup$ You're notation is a bit confusing, since you call a sphere a ball, and since $S^1$ and $S^2$ almost always denote the circle (in $R^2$) and the two-sphere (in $R^3$), respectively. $\endgroup$ – Hans Lundmark Mar 6 '11 at 15:59
14
$\begingroup$

Flip $n$ fair coins to pick an orthant—that is, to pick the signs of the coordinates of the point you are choosing. Now pick a point uniformly in the standard simplex, and flip the signs of its coordinates according to what the coins told you.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thanks. Step 1 was the trivial one. The reference for step 2 answers my question. $\endgroup$ – gappy Mar 7 '11 at 5:43
  • 7
    $\begingroup$ (+1) It is quite unfortunate that the authors of that paper have expended so much energy, and so recently(!), at promoting the Kraemer algorithm with $O(n \log n)$ runtime. A simpler $O(n)$ algorithm is as follows: Let $X_1, \ldots, X_n$ be iid unit exponential random variables and $S = \sum_i X_i$. Define $Y_i = X_i / S$. Then $(Y_1, \ldots, Y_n)$ is uniform in the standard simplex. Generating $X_i$ from standard uniforms is easy, just take $X_i = -\log(U_i)$ for $U_i \sim \mathcal U(0,1)$. For more general related issues, click here. $\endgroup$ – cardinal Sep 21 '11 at 19:18
9
$\begingroup$

Does the following work? (read the comment of joriki to convince yourself that the algorithm works --- thanks joriki)

Cut a line of length 1 into $n$ parts. Mathematically, throw $n-1$ random numbers between uniform between 0 and 1. Sort them to obtain $\mathbf{s}=(0,s_1, \dots, s_{n-1},1)$ with $s_i \leq s_j$; $\forall i <j$.

Take the point $\mathbf{x}= \left[\sigma_1(s_1- s_0), \sigma_2(s_2 - s_1), \dots, \sigma_n(s_n - s_{n-1})\right] \in S^1$. Here, $\sigma_i = \pm 1$ are randomly choosen. It is quite clear that $\mathbf{x}$ is on the unit sphere. The question remains: is it uniform on the sphere?

|cite|improve this answer
$\endgroup$
  • $\begingroup$ I'm reasonably convinced that it is uniform. Is it also simple enough? $\endgroup$ – Fabian Mar 6 '11 at 16:12
  • 1
    $\begingroup$ @Fabian: +1 -- This is brilliant! Can you offer any insight into how you came up with it? :-) I think uniformity follows by induction over $n$: The marginal distribution for $x_n$ is correct since $P(x_n>t)= (1-t)^{n-1}$, which is correct for a uniform distribution, and the correctness of the conditional distribution given $x_n$ follows from the induction hypothesis. $\endgroup$ – joriki Mar 6 '11 at 17:38
  • $\begingroup$ @joriki: I was just thinking how one can implement the "boundary condition" $\sum_i |x_i| =1$ without introducing bias in the individual variables (keeping the algorithm symmetric with respect to $i$). By the way, the sorting can of course be done on the fly (inserting the number at its appropriate place via bisection) and the only drawback of the algorithm is the memory usage of $n$ real numbers... $\endgroup$ – Fabian Mar 6 '11 at 17:53
  • $\begingroup$ @joriki: by the way it seems I don't get famous and this algorithm is known. The reference of Mariano Suárez-Alvarez calles it Krämers algorithm ;-) $\endgroup$ – Fabian Mar 6 '11 at 18:06
  • $\begingroup$ @Fabian: why would that be a problem? Your program will need to output a vector of n floats/doubles/... anyway, so you can't 'win' much in terms of memory allocation. $\endgroup$ – Gerben Mar 6 '11 at 18:27
1
$\begingroup$

Here's the method I proposed in the comment. Pick a random point $(x_1,\ldots,x_{n-1})$ from the $n-1$-dimensional hypercube. (This amounts to choosing $n-1$ real numbers uniformly in $[0,1]$)

Now, if $\sum_{i=1}^{n} x_i \geq 1$ throw the point out and start again, if not keep it.

Transform the point using the following affine map:

$$f:\mathbb{R}^{n-1}\to\mathbb{R}^{n}: (x_1,\ldots,x_{n-1}) \mapsto (x_1,\ldots,x_{n-1},1-\sum_{i=1}^{n} x_i) \; .$$

Now, similarly to what Fabian does, select $n$ signs $\sigma_i=\pm 1$ randomly and multiply each component of the vector you obtained by these signs.

As pointed out already by joriki and Gerben, for high dimensions, this method is very wasteful since a fraction $\frac{(n-1)!-1}{(n-1)!}$ points will have to be thrown out.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ As I already commented in response to your comment, this is rather wasteful, since it only uses one out of $(n-1)!$ points. Also I think it should read "$(n-1)$-dimensional hypercube". $\endgroup$ – joriki Mar 6 '11 at 16:27
  • $\begingroup$ Yeah, it's no good beyond dimension 3 or 4. $\endgroup$ – Raskolnikov Mar 6 '11 at 16:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.