how to generate random numbers that satisfy certain distribution

[Robert Kern]

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).

If you care about thread-safety, each thread should get its own Random instance 
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
  that is made terrible by our own mad attempt to interpret it as though it had
  an underlying truth."
   -- Umberto Eco