3
$\begingroup$

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!

$\endgroup$
6
  • 1
    $\begingroup$ Related: math.stackexchange.com/questions/67628/… $\endgroup$
    – cardinal
    Oct 2, 2011 at 13:39
  • $\begingroup$ @cardinal: Found the solution - look at my original question ;) $\endgroup$ Oct 2, 2011 at 14:06
  • $\begingroup$ 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. :) $\endgroup$
    – cardinal
    Oct 2, 2011 at 15:15
  • $\begingroup$ @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 :) $\endgroup$ Oct 2, 2011 at 15:36
  • 1
    $\begingroup$ 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. :) $\endgroup$
    – cardinal
    Oct 2, 2011 at 15:39

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$ Oct 2, 2011 at 12:35
  • 1
    $\begingroup$ @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. :) $\endgroup$
    – cardinal
    Oct 2, 2011 at 13:34
  • $\begingroup$ @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. ;) $\endgroup$ Oct 2, 2011 at 13:46
  • $\begingroup$ @Gortaur: I have a solution - please have a look at my question which is edited now. $\endgroup$ Oct 2, 2011 at 14:12

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .