By way of the central limit theorem, summing a dozen uniformly distributed variables (and normalizing) will give a a reasonable approximation to a normal distribution.

https://en.wikipedia.org/wiki/Central_limit_theorem

Very clever. My first instinct was that that would be inefficient, but addition is super cheap and the uniform distribution converges rathe quickly to normal. I like it a lot!

Even if LCG itself would be considered expensive, for quick prototyping needs this should be sufficient.

Once you have U([0,1]) you can generate anything you want. Just take Finv(u) where Finv is the inverse CDF of the thing you want. See [1] property 5.

https://en.wikipedia.org/wiki/Cumulative_distribution_functi...

Obviously. But as further reading exposes, the inverse is not always easily memorable -- or even easy to to compute! The inverse of the normal distribution, specifically, needs estimation of an integral to compute.

I guess it depends on what you assume are mathematical "primitives".

Log is an integral. Exp is an infinite series.

Normal cdf (erf) and its inverse (erfinv) are implemented almost everywhere exp and log are implemented. E.g. MKL and CUDA. Both can be computed via series or integral (like exp and log).

[1] https://en.wikipedia.org/wiki/Error_function

[2] https://people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf

Now that! Is much more constructive information. I'll look for those methods sooner in the future.

One can get two normally distributed random variables Z1 and Z2 from two uniformly distributed random variables U1 and U2, but the formulas are a little harder to remember:

https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform

There's good intuition for this. The key moral of the story is this:

> The normal distribution provides uniform sampling on the > unit n-sphere.

What does that mean? Assume one wants to sample a point uniformly from the unit circle. So we want to uniformly sample from the set { (x, y) : xx + yy = 1 }.

The naive method of uniformly sampling x over [0, 1] and then uniformly sampling y in [0, 1 - xx] doesn't work. It doesn't work because for a given change in dx, the change in dy is not uniform.

For example, if we are the point (x=0, y=1), in an infinitesimal region around this point, it looks like "flat land", so moving in dx does not change dy. That is, the tangent to the point (x=0, y=1) is horizontal.

On the other hand, if we are at the point (x=1, y=0), in an infinitesimal region around this point, it looks like a cliff face: all we have is vertical. So the dy change is "infinity". That is, the tangent to the point (x=1, y=0) is vertical.

Let's sample both x and y from the normal distribution. So the probability of getting x=r is e^{-rr}. The probability of getting y=s is e^{-ss}. The probability of getting some random point (x, y) = (r, s) is e^{-(rr+s*s)} = e^{-1} = constant.

So if we sample two points using the normal distribution, the points are uniformly distributed on the unit circle.

### Box–Muller transform

On the flip side, assume we have a uniformly distributed point on the circle. This means that the angle theta is uniformly distributed in [0, 2pi].

Given such a uniformly distributed point, we can generate two gaussian distributed numbers by taking their x and y component.

The intuition for why this works comes from the above:

- We saw that (x ~ normal, y ~ normal) => (theta ~ uniform) [as explained before].

- We can also show that (theta ~ uniform) => (x ~ normal, y ~ normal)

which is what the box-muller transform does.

Brilliant! I suspect this might help me remember the Box--Muller transform, which would be even better than relying on the CLT!