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).
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
More information about the Python-list
mailing list