Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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!

share|cite|improve this question
Related:… – cardinal Oct 2 '11 at 13:39
@cardinal: Found the solution - look at my original question ;) – Christian Ivicevic Oct 2 '11 at 14:06
Marsaglia's polar method is not a bad choice. A similar approach is the Box-Muller transform. It has the advantage that for every two uniform random variates input to the algorithm, one gets two independent standard normal variates out. So, for example, there is no need for the rejection sampling (i.e., loop) as used in the polar method. You "pay" for this extra convenience by needing to calculate $\sin$ and $\cos$, which are computationally a bit more expensive, but the price is cheap. :) – cardinal Oct 2 '11 at 15:15
@cardinal: I tried Box-Muller, however i was not able to shift the median like with Marsaglia! And if i think about the scenario: in the time i compute 1.000 uniform numbers i will need... 10 normal ones and i am confident, that i will be have with this :) – Christian Ivicevic Oct 2 '11 at 15:36
They're essentially the same method. If you get $(X_1,X_2)$ as the two outputs from each invocation of Box-Muller, then $a X_1 + b$ and $a X_2 + b$ are both $\mathcal N(b,a^2)$ and are independent of one another. At any rate, either method will work fine. :) – cardinal Oct 2 '11 at 15:39
up vote 2 down vote accepted

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.

share|cite|improve this answer

The normal distributed random variable can be just shifted in order to obtain the mean you need. E.g. in you algorithm it holds that $$ \sqrt{n}(S_n-\mu)+a \to \mathcal N(a,\sigma^2). $$

But if you just want to obtain the normal distributed random variable from the uniform distribution, why don't just use a random variable $U$ which is uniformly distributed on $[0,1]$ - then $F^{-1}(U)$ will be distributed as a random variable with cumultaive distribution function $F$.

For example, for normal random variable, $$ F(x;\mu,\sigma^2) = \frac1{\sigma\sqrt{2\pi}}\int\limits_{-\infty}^x\mathrm e^{-(t-\mu)^2/2\sigma^2}\mathrm dt $$ and I believe this function you can easily calculate numerically as well as its inverse. Moreover, if you are interested in discrete value, you may want to use Binomial Distribution which is very similar to the normal.

share|cite|improve this answer
I must confess, that adding a scalar $a$ was my first idea, however my algorithm seems buggy and it doesn't work like expected. Furthermore, debugging my code revealed that i am getting sometimes numbers $\leq 0$ which is weird. Any suggestions what i can do? – Christian Ivicevic Oct 2 '11 at 12:35
@Christian: Of course you're getting negative numbers, and about half of the time! The distribution of $S - \mu$ is symmetric (i.e., even) about zero. In fact, it's a symmetric function of piecewise polynomials. By the way, this is an incredibly inefficient way to simulate an approximately normal random variable. :) – cardinal Oct 2 '11 at 13:34
@cardinal: I need this only 10 times at the beginning of the game and maybe in really special circumstances and therefore efficiency is not that important. ;) – Christian Ivicevic Oct 2 '11 at 13:46
@Gortaur: I have a solution - please have a look at my question which is edited now. – Christian Ivicevic Oct 2 '11 at 14:12

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.