# how to generate random numbers that satisfy certain distribution

Robert Kern robert.kern at gmail.com
Mon Jan 25 02:33:42 CET 2010

```On 2010-01-23 21:30 , Steven D'Aprano wrote:
> On Sat, 23 Jan 2010 14:10:10 -0800, Paul Rubin wrote:
>
>> Steven D'Aprano<steve at REMOVE-THIS-cybersource.com.au>  writes:
>>> The Box-Muller transform is reasonably simple, but you can't be serious
>>> that it is simpler than adding twelve random numbers and subtracting
>>> six!
>>
>> If you want a normal distribution, using the Box-Muller transform is
>> simpler because it spares you the complication of figuring out whether
>> the 12-trial binomial approximation is close enough to produce reliable
>> results for your specific application, which you obviously have to do if
>> you are using the approximation for anything serious.
>
> But using Box-Miller gives you the complication of figuring out whether
> you care about being thread safety, which you have to do if you're doing
> anything serious. (See the comments in the random module for the gauss
> method).

anyways.

> By the way, does anyone know why there is no Poisson random numbers in
> the module? The implementation is quite simple (but not as simple as the
> Linux kernel *wink*):
>
>
> def poisson(lambda_=1):
>      L = math.exp(-lambda_)
>      k = 0
>      p = 1
>      while 1:
>          k += 1
>          p *= random.random()
>          if p<= L:
>              break
>      return k-1
>
>
> although this can be improved upon for large values of lambda_.

Where large is about 10. That's where my implementation in numpy.random switches
over to a more complicated, but more efficient implementation. It does require a
log-gamma special function implementation, though.

#define LS2PI 0.91893853320467267
#define TWELFTH 0.083333333333333333333333
long rk_poisson_ptrs(rk_state *state, double lam)
{
long k;
double U, V, slam, loglam, a, b, invalpha, vr, us;

slam = sqrt(lam);
loglam = log(lam);
b = 0.931 + 2.53*slam;
a = -0.059 + 0.02483*b;
invalpha = 1.1239 + 1.1328/(b-3.4);
vr = 0.9277 - 3.6224/(b-2);

while (1)
{
U = rk_double(state) - 0.5;
V = rk_double(state);
us = 0.5 - fabs(U);
k = (long)floor((2*a/us + b)*U + lam + 0.43);
if ((us >= 0.07) && (V <= vr))
{
return k;
}
if ((k < 0) ||
((us < 0.013) && (V > us)))
{
continue;
}
if ((log(V) + log(invalpha) - log(a/(us*us)+b)) <=
(-lam + k*loglam - loggam(k+1)))
{
return k;
}

}

}

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma