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\}$?

  • 1
    Have you tried the same procedure but normalizing using the 1-norm?2011-03-06
  • 1
    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.$2011-03-06
  • 3
    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$?2011-03-06
  • 0
    [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.]2011-03-06
  • 4
    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.2011-03-06
  • 0
    How about the following: the "1-ball" in $n$ dimensions is consisting of $2^n$ "flat $n-1$-faces" in $2^n$ $2^n$-ants. Pick uniformly distributed points in a $n-1$ dimensional hypercube. Throw out the points with sum of coordinates bigger than 1. For each of these points also randomly select $n$-vectors consisting of $1$ or $-1$ as coordinates. Now, by an affine transformation, map the points you selected to points on the face in the "positive" $2^n$-ant. Multiply now the points with their corresponding sign-vector coordinate-wise. Presto! (You can tell I'm not a programmer.)2011-03-06
  • 0
    @Raskolnikov: the problem resides with 'throwing out points'. As you increase the dimension of your space, you quickly throw out almost all points, making it unfeasible to practically do this.2011-03-06
  • 0
    @Raskolnikov: This only uses $1/n!$ of the points -- that's rather wasteful :-)2011-03-06
  • 0
    @Gerben, @joriki: Told you I'm not a good programmer. ;p2011-03-06
  • 0
    @joriki correct. @lundmark correct. Sorry for the sloppy language. It was edited. @gerben correct. I thought of that procedure, but it's highly inefficient.2011-03-07

3 Answers 3

14

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.

  • 0
    Thanks. Step 1 was the trivial one. The reference for step 2 answers my question.2011-03-07
  • 7
    (+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](http://en.wikipedia.org/wiki/Dirichlet_distribution).2011-09-21
8

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

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?

  • 0
    I'm reasonably convinced that it is uniform. Is it also simple enough?2011-03-06
  • 1
    @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.2011-03-06
  • 0
    @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...2011-03-06
  • 0
    @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 ;-)2011-03-06
  • 0
    @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.2011-03-06
  • 0
    @Gerben: oh, that is true... So no drawback...2011-03-06
  • 0
    @Fabian: (+1) This algorithm is quite nice and clever. One drawback is that it is $O(n\log n)$, whereas there are simple $O(n)$ approaches. See my comment to Mariano.2011-09-21
1

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.

  • 0
    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".2011-03-06
  • 0
    Yeah, it's no good beyond dimension 3 or 4.2011-03-06