# Which distributions can you sample if you can sample a Gaussian?

Explicitly: You have a computer that is able to pick a real number at random according to the normal distribution: $\mathcal{N}(0,1) = \frac{1}{\sqrt{2\pi}}e^{-x^2/2}$. Which distributions can this computer sample from, provided that your program must terminate after finitely many steps (not almost surely, but logically must terminate)?

For a simpler example, if your computer can sample a 'fair coin', heads 1/2 and tails 1/2, then it can also sample any binary event with probability a dyadic rational: $p/2^n$, $0 \leq p \leq 2^n$ an integer. It can't sample an event with probability $1/3$. The program "Flip two coins, $HH$ is a success, and on $TT$ repeat the experiment" is not allowed because it is not guaranteed to terminate (even though it will terminate with probability $1$).

Since it might be hard to characterize 'all possible' distributions, here's a concrete question: Can you sample from the Haar measure on any compact matrix Lie group?

EDIT: The idea of a random real number has more philosophical and practical subtleties than I realized when posting this question. I encourage commenters to answer based on varying interpretations, but to get the ball rolling how about these restrictions:

1) Given random real numbers x and y, the computer can calculate x + y, xy, and x/y (provided y =/= 0) in one step. (You can also multiply/add/divide by rational numbers)

2) Given a real random number x, the computer can decide in one step if x > 0, x = 0, or x < 0 and act accordingly.

I know this might not accurately model how computers or computer models work, so like I said feel free to play with these assumptions.

• Can you explain the sense in which the computer has the real number? What is your computational formalism here? Should we imagine that there are some registers (as many as desired?) that hold the "random" real numbers? Can we, for example, add them together "in finitely many steps"? How are the algorithms manipulating real numbers? Are you using a Turing machine model, the BSS model, or what? – Joel David Hamkins Apr 30 '14 at 19:45
• I edited with some assumptions. After wiki-ing, it looks like I'm effectively using the BSS model. – Alex Zorn Apr 30 '14 at 20:37
• If you add in the constant $\pi$ and the ability to perform elementary functions, then being able to sample from the normal distribution is equivalent to being able to sample from a uniform distribution because the Box-Mueller transform only requires those (in both directions). en.wikipedia.org/wiki/Box%E2%80%93Muller_transform – Douglas Zare Apr 30 '14 at 23:21
• You shouldn't use the BSS model, it has unrealistic assumptions, such as that two real numbers can be compared for equality. – Andrej Bauer May 1 '14 at 7:49
• @AlexZorn Isn't it more natural to require only that the algorithms almost surely terminate, rather than logically always? Your requirement goes against what I see as a fundamental philosophy for this kind of context, where we should ignore sets of measure zero. Requiring the algorithms to work only almost surely avoids many irritating issues and will enable you to sample additional distributions, as you point out already in the question. – Joel David Hamkins May 1 '14 at 13:32

This is by no means a complete classification, although I'm surprised that no-one has mentioned that we can easily construct the most well-known continuous distributions from rational functions of $$N(0,1)$$ variables.

Let $$\{ X_1, X_2, X_3, \dots \}$$ be independent $$N(0,1)$$ random variables. Then we have:

• $$\textrm{Gamma}(\frac{n}{2},\frac{1}{2}) = X_1^2 + \dots + X_n^2$$
• $$\textrm{Beta}(\frac{n}{2},\frac{m}{2}) = \dfrac{X_1^2 + \dots + X_n^2}{X_1^2 + \dots + X_{n+m}^2}$$
• $$\textrm{Uniform}[0,1] = \dfrac{X_1^2 + X_2^2}{X_1^2 + X_2^2 + X_3^2 + X_4^2}$$ (special case of Beta)

You can also get the Student's $$t$$-distribution, amongst other things (slightly more complicated, but still a rational function of independent $$N(0,1)$$ random variables).

Anyway, by adding and multiplying by constants, we can rescale those distributions:

• $$\textrm{Uniform}[a,b] = \dfrac{a(X_1^2 + X_2^2) + b(X_3^2 + X_4^2)}{X_1^2 + X_2^2 + X_3^2 + X_4^2}$$
• $$\textrm{Gamma}(\frac{n}{2},\lambda) = 2\lambda(X_1^2 + \dots + X_n^2)$$
• $$\textrm{Exponential}(\lambda) = 2\lambda(X_1^2 + X_2^2)$$ (special case of Gamma)

I think that being able to construct the uniform distribution is probably my favourite of these results, not least because it's not immediately obvious from the construction that it should be uniform.

EDIT: On deeper investigation, it transpires that the construction of the uniform distribution in terms of four Gaussians has an elegant interpretation in the context of quantum computation.

• This is a great answer. Any intuition for why you need 4 squared Gaussians for the uniform? – Bjørn Kjos-Hanssen May 3 '14 at 5:03
• Yes, we can think of it as choosing a random quaternion q according to a spherically-symmetric (Gaussian) distribution. Then the formula for Uniform[-1, 1] is just one coordinate of the image of q/|q| under the Hopf map, and we're done by Archimedes' theorem that the orthogonal projection of a uniform distribution over S^2 is a U[-1, 1] distribution. – Adam P. Goucher Jul 5 '14 at 18:39

Just want to observe that if $X$ is $\mathcal N(0,1)$ and $Y$ has any distribution with a strictly increasing* cdf $F_Y$ such that $F_Y^{-1} \circ F_X$ is computable (in whatever model of computability you are working with), then each one is sampleable from the other.

Indeed suppose we have a $\mathcal N(0,1)$ random variable $X$ with cdf $F_X$.

Note the random variable $F_Y(Y)$ is uniform on $[0,1]$.

Suppose we obtain a sample value $x$ from $X$. Then $F_Y^{-1}(F_X(x))$ is a sample value from $Y$: $$\mathbb P(Y\le F_Y^{-1}(F_X(x))) = \mathbb P(F_Y(Y)\le F_X(x))=F_X(x)=\mathbb P(X\le x).$$ *Actually it's enough that $Y$ have a strictly increasing cdf on an interval that it a.s. belongs to.

However, $x\mapsto e^x$ is not BSS-computable, so there may not be that many such $Y$'s if you use that model.

Lie groups are made up of spheres and it is easy to sample from spheres, so it is easy to sample from Haar measure. We can sample from $S^{n-1}\subset\mathbb R^n$ by taking $n$ normal samples, assuming that they aren't all zero, and normalizing the vector to lie on the unit sphere, using square roots. $SO(n)$ acts on $S^{n-1}$, preserving the standard measure, with stabilizer conjugate to $SO(n-1)$. By induction we can sample from $SO(n-1)$, so we can sample from $SO(n)$. Other compact connected groups have similar structures. If we want to sample from disconnected groups, we must divide our sample space into $n$ regions of equal measure. This may be tricky for the difficult distribution, but it is easy for the circle, given trig functions.