Transform Uniform Distribution to Normal Using Lindeberg-Levy CLT – Probability

limitsprobabilityprobability distributionsprobability theory

Currently i am developing a game which is based on many computations of random values and therefore i have implemented many algorithms like the Mersenne-Twister etc. Unfortunately, all generators return uniformly distributed values and i want to modify this to make the min/max (and their surrounding) values rarer than the other ones. Researches brought me up to the point that i have been using now the Lindeberg–Lévy Central Limit Theorem:

Theorem. Suppose {$X_i$} is a sequence of iid random variables with $\mathbb{E}[X_i] = \mu$ and $Var[X_i] = \sigma^2$. Let $S_n=\frac{1}{n}(X_1+X_2+\ldots+X_n)$. Then as $n$ approaches infinity, the random variable $\sqrt{n}(S_n − \mu)$ converges in distribution to a normal $\mathcal{N}(0, \sigma^2)$:
$$\sqrt{n}\cdot\left(\left(\frac{1}{n}\sum\limits_{i=1}^nX_i\right)-\mu\right)\overset{d}{\rightarrow}\mathcal{N}(0, \sigma^2)$$

My C++ code works fine… however the median is set to 0 and i want to know how i can change the median and the range in combination to get e.g. values from $[1;100]$ as a result will have the property, that the median ones (50 and surrounding) will have a higher probability to be returned than [1;5] or [96;100]. (Hope you understand what i am thinking right now!)

Concrete example

Throwing 10.000 "unfair" dice currently leads to this distribution: (8836,729,341,86,8,1) however i want something that looks like this: (18,341,4418,4568,341,56)

Current algorithm (old one based on the theorem!)

    n <- 32;
    S <- 0.0;
    for i = 0 to n
        S <- S + mersenneTwister(min, max)
    S <- S / n
    µ <- (min + max) / 2.0
    randval <- sqrt(n) * (S - µ)
    return max(min, floor(randval) mod max)

It seems to me, that something it not implemented correctly…

Solution (i am not yet able to answer my own questions…)

Finally i have found a very simple solution which suits my needs and is no overkill at all.

Marsaglia polar method (more information here)

  1. Generate two iid random values $u_1$ and $u_2$ in [-1;1].
  2. Compute $q = u_1^2 + u_2^2$. In case that $q=0$ or $q > 1$ return to step 1.
  3. Compute $p=\sqrt{-2\cdot \ln(q) / q}$.
  4. $x_{1/2} = u_{1/2}\cdot p$ represent two independant normal distributed random numbers.
  5. If $x\sim\mathcal{N}(0,1)$, then $a\cdot x+b\sim\mathcal{N}(b,a^2)$.

Exactly what i need! For throwing my dice i let $a=1$ and $b=\mu=(min+max)/2$.

@cardinal: Thanks for this idea!

Best Answer

The Marsaglia Polar Method is an algorithmic improvement on the Box-Muller Transform. The Box-Muller Transform uses two uniform random variables, $(u,v)$, in $[0,1]\times[0,1]$ to generate two normally distributed variables in $\mathbb{R}\times\mathbb{R}$: $$ (\cos(2\pi u),\sin(2\pi u))\sqrt{-2\log(v)} $$ The Marsaglia Polar Method uses the fact that for a uniformly distributed point, $(x,y)$, in the unit disk, $x^2+y^2$ is uniformly distributed in $[0,1]$. Thus, $x^2+y^2$ can take the place of $v$ in the Box-Muller Transform. Since $\frac{(x,y)}{\sqrt{x^2+y^2}}$ is uniformly distributed on the unit circle, it can take the place of $(\cos(2\pi u),\sin(2\pi u))$ in the Box-Muller Transform. With these substitutions, we get that $$ \frac{(x,y)}{\sqrt{x^2+y^2}}\sqrt{-2\log(x^2+y^2)} $$ is a normally distributed random variable. The Marsaglia Polar Method uses a Monte-Carlo method to generate a uniformly distributed $(x,y)$ in the unit disk; that is, it picks points in $[-1,1]\times[-1,1]$ until one lands in the unit disk (which happens with probability $\frac{\pi}{4}$). To prevent illegal values, the unlikely value $(0,0)$ is ignored.

I say that the Marsaglia Polar Method is an algorithmic improvement on the Box-Muller Transform because it avoids the computation of a trigonometric function or two, but at the cost of generating $\frac{4}{\pi}$ uniform pairs, on average, per each normal pair.

Related Question