jdscheff/basketball-gm — pull request #14 by kgilbert-cmu (2013-07-30 10:38:02+00:00)
Title: Random Gaussians are decidedly not Gaussian

[This](http://i.imgur.com/esA4Qhn.jpg) is an estimate of the quantile-comparison for your `random.gauss()` method compared to an actual random Gaussian distribution (comparison was made against the library standard `rnorm` function in the R Project). You'll notice that the far edges of the line curves away from the straight `y=x` line. That's not supposed to happen. This indicates that in the tail of the distribution, your algorithm is producing skewed results. I simulated about 100,000 draws from each, and your current method to select Gaussians has **extremely** short tails. I think it cuts about an entire standard deviation off. After that many replicates, I should have seen about 270 random draws that had an absolute value greater than 3.0 (so more than 3 standard deviations away). There were none. The closest was 2.97.

And in fact, I know why.

The website you linked in the comment directed the reader to some code that reproduces the [Uniform Sum Distribution](http://mathworld.wolfram.com/UniformSumDistribution.html). The concept is that if you take N uniform random variables and add them together, then the distribution of the sum of those uniforms is actually a Gaussian, with a mean and standard deviation and everything.

Unfortunately, it only *approaches* a Gaussian, and even then only if you add a whole bunch together. The current usage on `jdscheff:master` takes the sum of merely 2 random uniforms, and that's not nearly enough. Then, to really mess with things, it's like you added an additional `sigma` worth of the uniforms, but with none of the entropy to show for it.

I've replaced the entire function with an implementation of Marsaglia's polar method. Draw a square (Cartesian coordinate plane) and inscribe the unit circle on it. Pick a random point inside the square (call it (X,Y) or (Z_1, Z_2)) and only use it if it fits inside the bounds of the circle. Find the distance between the center of the circle and this point, call that S. Then, for both Z, the number `Z * sqrt(-2 * log(S) / S)` is a random Gaussian centered at 0 with standard deviation 1. We then scale by the target standard deviation and add to target mean to get desired Gaussian distribution.

I simulated an additional 100,000 random draws from the new algorithm. [The quantiles plot looks like this](http://i.imgur.com/WOBf8uK.jpg)

Ah, much better.