Re: random numbers according to user defined distribution ??

2008-08-08 Thread Robert Kern

sturlamolden wrote:

Alex wrote:


I wonder if it is possible in python to produce random numbers
according to a user defined distribution?
Unfortunately the random module does not contain the distribution I
need :-(


There exist some very general algorithms to generate random numbers
from arbitrary distributions.

The most notable of these are Markov Chain Monte Carlo, e.g. the
Metropolis-Hastings algorithm. It is very easy to implement in any
programming language. The nature MCMC algorithms makes it inefficient
when implemented in pure Python. But you can get tremendous speedup by
simulating multiple Markov chains in parallel, by means of vectorizing
with NumPy.

A relative of Metropolis-Hastings which may also be applicable to your
problem is pure rejection sampling. It is far less efficient, but
produces no autocorrelation in the samples.


I don't know. I think the certainty that rejection sampling actually gives you 
the desired distribution as opposed to MH's uncertainty is very much a 
worthwhile tradeoff for univariate distributions. Sure, you throw away fewer 
samples, but you know that MH isn't throwing away samples that it ought to. 
Personally, I view MH as a last resort when the dimensionality gets too large to 
do anything else. But then, that's just my opinion.


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

--
http://mail.python.org/mailman/listinfo/python-list


Re: random numbers according to user defined distribution ??

2008-08-08 Thread Jeremy Sanders
Alex wrote:

 I wonder if it is possible in python to produce random numbers
 according to a user defined distribution?
 Unfortunately the random module does not contain the distribution I
 need :-(

Have you looked at the numpy random number module? It seems to have quite a
lot of distributions. See help(numpy.random).

Jeremy

-- 
Jeremy Sanders
http://www.jeremysanders.net/
--
http://mail.python.org/mailman/listinfo/python-list


Re: random numbers according to user defined distribution ??

2008-08-07 Thread Raymond Hettinger
On Aug 6, 3:02 pm, Alex [EMAIL PROTECTED] wrote:
 I wonder if it is possible in python to produce random numbers
 according to a user defined distribution?
 Unfortunately the random module does not contain the distribution I
 need :-(

Sure there's a way but it won't be very efficient.  Starting with an
arbitrary probability density function over some range, you can run it
through a quadrature routine to create a cumulative density function
over that range.  Use random.random() to create a uniform variate x.
Then use a bisecting search to find x in the cumulative density
function over the given range.

from __future__ import division
from random import random

def integrate(f, lo, hi, steps=1000):
dx = (hi - lo) / steps
lo += dx / 2
return sum(f(i*dx + lo) * dx for i in range(steps))

def make_cdf(f, lo, hi, steps=1000):
total_area = integrate(f, lo, hi, steps)
def cdf(x):
assert lo = x = hi
return integrate(f, lo, x, steps) / total_area
return cdf

def bisect(target, f, lo, hi, n=20):
'Find x between lo and hi where f(x)=target'
for i in range(n):
mid = (hi + lo) / 2.0
if target  f(mid):
hi = mid
else:
lo = mid
return (hi + lo) / 2.0

def make_user_distribution(f, lo, hi, steps=1000, n=20):
cdf = make_cdf(f, lo, hi, steps)
return lambda: bisect(random(), cdf, lo, hi, n)

if __name__ == '__main__':
def linear(x):
return 3 * x - 6
lo, hi = 2, 10
r = make_user_distribution(linear, lo, hi)
for i in range(20):
print r()

--
Raymond
--
http://mail.python.org/mailman/listinfo/python-list


Re: random numbers according to user defined distribution ??

2008-08-07 Thread Alex
Thanks for the many answers.

So basically I have to get the inverse of the CDF and use this to
transform my uniformly distributed random numbers. If my desired
distribution is simple I can get an analytical solution for the
inverse, otherwise I have to use numerical methods.

Okay, things are now much clearer.

Many thanks,

 Alex
--
http://mail.python.org/mailman/listinfo/python-list


Re: random numbers according to user defined distribution ??

2008-08-07 Thread Robert Kern

Alex wrote:

Thanks for the many answers.

So basically I have to get the inverse of the CDF and use this to
transform my uniformly distributed random numbers. If my desired
distribution is simple I can get an analytical solution for the
inverse, otherwise I have to use numerical methods.

Okay, things are now much clearer.


It's worth noting that unless if the PPF (the inverse of the CDF) is very 
straightforward, this method is not very good. The numerical errors involved 
cause very poor results in the tails of many distributions. Typically, rejection 
sampling, if done well, will work much better. There are techniques for doing 
this on nearly-arbitrary distributions:


  http://statmath.wu-wien.ac.at/projects/arvag/index.html

If you implement any of these techniques in Python, I would love to see them.

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

--
http://mail.python.org/mailman/listinfo/python-list


Re: random numbers according to user defined distribution ??

2008-08-07 Thread sturlamolden
Alex wrote:

 I wonder if it is possible in python to produce random numbers
 according to a user defined distribution?
 Unfortunately the random module does not contain the distribution I
 need :-(

There exist some very general algorithms to generate random numbers
from arbitrary distributions.

The most notable of these are Markov Chain Monte Carlo, e.g. the
Metropolis-Hastings algorithm. It is very easy to implement in any
programming language. The nature MCMC algorithms makes it inefficient
when implemented in pure Python. But you can get tremendous speedup by
simulating multiple Markov chains in parallel, by means of vectorizing
with NumPy.

A relative of Metropolis-Hastings which may also be applicable to your
problem is pure rejection sampling. It is far less efficient, but
produces no autocorrelation in the samples.

If you can generate random deviates from the marginal distributions,
but need to sample the joint distribution, look into using the Gibbs
sampler. It is an efficient version of Metropolis-Hastings for this
special case.

http://en.wikipedia.org/wiki/Metropolis-Hastings_algorithm
http://en.wikipedia.org/wiki/Rejection_sampling
http://en.wikipedia.org/wiki/Gibbs_sampling

--
http://mail.python.org/mailman/listinfo/python-list


random numbers according to user defined distribution ??

2008-08-06 Thread Alex
Hi everybody,

I wonder if it is possible in python to produce random numbers
according to a user defined distribution?
Unfortunately the random module does not contain the distribution I
need :-(

Many thanks

 axel
--
http://mail.python.org/mailman/listinfo/python-list


Re: random numbers according to user defined distribution ??

2008-08-06 Thread Dominic van Berkel
On Thursday 07 August 2008 00:02, Alex [EMAIL PROTECTED] wrote:

 Hi everybody,
 
 I wonder if it is possible in python to produce random numbers
 according to a user defined distribution?
 Unfortunately the random module does not contain the distribution I
 need :-(
 
 Many thanks
 
  axel
I'm not aware of any module with that specific function, but it's
algorithmically not too complex I'd think. If you're writing an application
that does this I'll assume that you have a basic gist of how to implement
it ;). It's been a while since I messed with the subject, so I'm not
getting any further than graphs in my head right now. Good luck!

-- 
Dominic van Berkel
Bi-la Kaifa
--
http://mail.python.org/mailman/listinfo/python-list


Re: random numbers according to user defined distribution ??

2008-08-06 Thread Steven D'Aprano
On Wed, 06 Aug 2008 15:02:37 -0700, Alex wrote:

 Hi everybody,
 
 I wonder if it is possible in python to produce random numbers according
 to a user defined distribution? Unfortunately the random module does not
 contain the distribution I need :-(


This is a strange question. Of course you can -- just write a function to 
do so! Here's some easy ones to get you started:

from __future__ import division
import random, maths

def unbounded_rand(p=0.5):
Return a random integer between 0 and infinity.
if not (0  p = 1):
raise ValueError
n = 0
while random.random()  p:
n += 1
return n

def pseudonorm():
Return a random float with a pseudo-normal distribution.

The probability distribution is centered at 0 and bounded 
by -1 and +1.

return (sum([random.random() for i in range(6)])-3)/3

def triangular(min=0, max=1, mode=0.5):
Return a random float in the range (min, max) inclusive
with a triangular histogram, and the peak at mode.

u = random.random()
if u = (mode-min)/(max-min):
return min + math.sqrt(u*(max-min)*(mode-min))
else:
return max - math.sqrt((1-u)*(max-min)*(max-mode))

def linear():
Return a random float with probability density 
function pdf(x)=2x.

return math.sqrt(random.random())



There's no general way to create a random function for an arbitrary 
distribution. I don't think there's a general way to *describe* an 
arbitrary random distribution. However, there are some mathematical 
techniques you can use to generate many different distributions. Google 
on transformation method and rejection method.

If you have a specific distribution you are interested in, and you need 
some help, please ask.



-- 
Steven
--
http://mail.python.org/mailman/listinfo/python-list


Re: random numbers according to user defined distribution ??

2008-08-06 Thread Dan Bishop
On Aug 6, 8:26 pm, Steven D'Aprano [EMAIL PROTECTED]
cybersource.com.au wrote:
 On Wed, 06 Aug 2008 15:02:37 -0700, Alex wrote:
  Hi everybody,

  I wonder if it is possible in python to produce random numbers according
  to a user defined distribution? Unfortunately the random module does not
  contain the distribution I need :-(

 This is a strange question. Of course you can -- just write a function to
 do so! Here's some easy ones to get you started:

...

 There's no general way to create a random function for an arbitrary
 distribution. I don't think there's a general way to *describe* an
 arbitrary random distribution.

What about the quantile function?
--
http://mail.python.org/mailman/listinfo/python-list


Re: random numbers according to user defined distribution ??

2008-08-06 Thread Paul Rubin
Alex [EMAIL PROTECTED] writes:
 I wonder if it is possible in python to produce random numbers
 according to a user defined distribution?

That can mean a lot of things.  The book Numerical Recipes (there
are editions for various languages, unfortunately not including Python
last time I looked) has some discussion about how to do it.

--
http://mail.python.org/mailman/listinfo/python-list


Re: random numbers according to user defined distribution ??

2008-08-06 Thread Steven D'Aprano
On Wed, 06 Aug 2008 21:09:30 -0700, Dan Bishop wrote:

 There's no general way to create a random function for an arbitrary
 distribution. I don't think there's a general way to *describe* an
 arbitrary random distribution.
 
 What about the quantile function?


Well, sure, if you can write down the quantile function, c.d.f or p.d.f. 
of a distribution, I suppose that counts as describing it, in some sense. 
But even if we limit ourselves to distributions which are actually 
useful, as opposed to arbitrary distributions that can't be described in 
terms of any known mathematical function, there are serious practical 
difficulties. I quote from the Wikipedia article on quantile functions:

The quantile functions of even the common distributions are relatively 
poorly understood beyond the use of simple lookup tables, which is at 
odds with their importance in Monte Carlo sampling, where a sample from a 
given distribution may be obtained in principle by applying its quantile 
function to a sample from a uniform distribution. The exponential case 
above is one of the very few distributions where there is a simple 
formula.

http://en.wikipedia.org/wiki/Quantile_function


-- 
Steven
--
http://mail.python.org/mailman/listinfo/python-list