Prev: how can i know if a python object have a attribute such as 'attr1'?
Next: ISO module for binomial coefficients, etc.
From: Paul Rubin on 23 Jan 2010 15:26
thinke365 <thinke365(a)gmail.com> writes:
> of course i have tried random package, but can this package generate random
> sequence that satisfy possion distribution , normal distribution and uniform
> distribution

Did you look at the documentation?
From: Paul Rubin on 23 Jan 2010 15:29
Peter Chant <peteRE(a)MpeteOzilla.Vco.ukE> writes:
> I remeber being told that adding up 12 random numbers in the range 0-1
> (which is what most computer random number genertors at the time chucked
> out) and subtracted 6 gives a pretty good normal distribution. I think I
> did try it once and it failed, but I must have done something odd.

That gives you a binomial distribution on 12 trials, which approximates
a normal distribution when the number of trials is large. 12 isn't too
bad. But there's a simpler way, the Box-Muller transform, that gives
you a pair drawn from a precisely normal distribution from two uniform
random samples:

http://en.wikipedia.org/wiki/Box-Muller_transform
From: Paul Rubin on 23 Jan 2010 17:10
Steven D'Aprano <steve(a)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. It also involves
writing less code than that list comprehension, since it is already
implemented in the random module so you can just call it directly.
From: Robert Kern on 24 Jan 2010 20:33
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(a)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

First  |  Prev  | 
Pages: 1 2
Prev: how can i know if a python object have a attribute such as 'attr1'?
Next: ISO module for binomial coefficients, etc.