random numbers according to user defined distribution ??

Discussion in 'Python' started by Alex, Aug 6, 2008.

  1. Alex

    Alex Guest

    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
     
    Alex, Aug 6, 2008
    #1
    1. Advertising

  2. On Thursday 07 August 2008 00:02, 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 :-(
    >
    > 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"
     
    Dominic van Berkel, Aug 6, 2008
    #2
    1. Advertising

  3. 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
     
    Steven D'Aprano, Aug 7, 2008
    #3
  4. Alex

    Dan Bishop Guest

    On Aug 6, 8:26 pm, Steven D'Aprano <st...@REMOVE-THIS-
    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?
     
    Dan Bishop, Aug 7, 2008
    #4
  5. Alex

    Paul Rubin Guest

    Alex <> 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.
     
    Paul Rubin, Aug 7, 2008
    #5
  6. 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
     
    Steven D'Aprano, Aug 7, 2008
    #6
  7. On Aug 6, 3:02 pm, 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 :-(


    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
     
    Raymond Hettinger, Aug 7, 2008
    #7
  8. Alex

    Alex Guest

    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
     
    Alex, Aug 7, 2008
    #8
  9. Alex

    Robert Kern Guest

    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
     
    Robert Kern, Aug 7, 2008
    #9
  10. Alex

    sturlamolden Guest

    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
     
    sturlamolden, Aug 8, 2008
    #10
  11. Alex

    Robert Kern Guest

    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
     
    Robert Kern, Aug 8, 2008
    #11
  12. 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/
     
    Jeremy Sanders, Aug 8, 2008
    #12
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Oodini
    Replies:
    1
    Views:
    1,796
    Keith Thompson
    Sep 27, 2005
  2. globalrev
    Replies:
    4
    Views:
    782
    Gabriel Genellina
    Apr 20, 2008
  3. petertwocakes

    Liinear distribution of random numbers?

    petertwocakes, Jan 17, 2009, in forum: C Programming
    Replies:
    32
    Views:
    892
    CBFalconer
    Jan 25, 2009
  4. thinke365
    Replies:
    11
    Views:
    2,295
    Robert Kern
    Jan 25, 2010
  5. VK
    Replies:
    15
    Views:
    1,217
    Dr J R Stockton
    May 2, 2010
Loading...

Share This Page