How to generate geometric random numbers?

Discussion in 'Python' started by MyInfoStation, Jul 23, 2006.

  1. Hi all,

    I am a newbie to Python and would like to genereate some numbers
    according to geometric distribution. However, the Python Random package
    seems do not have implemented functionality. I am wondering is there
    exist any other libraries that can do this job?

    Thanks a lot,

    Da
     
    MyInfoStation, Jul 23, 2006
    #1
    1. Advertisements

  2. The usual way is to generate standard random numbers (linear distribution)
    and then apply whatever transformation you need to generate the desired
    distribution.

    Gerhard
     
    Gerhard Fiedler, Jul 23, 2006
    #2
    1. Advertisements

  3. MyInfoStation

    Robert Kern Guest

    That only works if there is such a transformation.

    The geometric distribution and many others have been implemented in numpy:

    http://www.scipy.org/NumPy

    In [1]: from numpy import random

    In [2]: random.geometric(0.5, size=100)
    Out[2]:
    array([1, 5, 2, 3, 1, 1, 2, 3, 1, 1, 2, 1, 1, 1, 2, 2, 3, 3, 1, 1, 1, 1, 1,
    2, 2, 1, 1, 1, 1, 1, 3, 1, 1, 3, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 3, 1,
    4, 1, 1, 1, 2, 1, 2, 3, 2, 1, 1, 1, 1, 1, 3, 1, 1, 2, 6, 1, 1, 3, 2,
    1, 1, 2, 1, 1, 7, 2, 1, 1, 2, 1, 1, 2, 4, 1, 2, 1, 4, 2, 1, 1, 2, 1,
    4, 2, 1, 1, 3, 1, 3, 1])

    --
    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, Jul 23, 2006
    #3
  4. Thanks a lot. I will try it out.

    But I am still surprised because the default Random package in Python
    can generate so few discrete random distritbuions, while it can
    generate quite a few continuous distribution, including some not very
    common one.

    Da
     
    MyInfoStation, Jul 24, 2006
    #4
  5. MyInfoStation

    Paul Rubin Guest

    It looks pretty simple to transform the uniform distribution to the
    geometric distribution. The formula for its cdf is pretty simple:

    cdf(p,n) = (1-p)**(n-1)*p

    For fixed p, if the cdf is c, we get (unless I made an error),

    n = log(c, 1-p) - 1

    So choose a uniform point c in the unit interval, run it through that
    formula, and round up to the nearest integer.

    See http://en.wikipedia.org/wiki/Geometric_distribution
    for more about the distribution.
     
    Paul Rubin, Jul 24, 2006
    #5
  6. MyInfoStation

    Paul Rubin Guest

    I meant n = log(c/p, 1-p) - 1
    sorry.
     
    Paul Rubin, Jul 24, 2006
    #6
  7. MyInfoStation

    Robert Kern Guest

    import random
    from math import ceil, log

    def geometric(p):
    """ Geometric distribution per Devroye, Luc. _Non-Uniform Random Variate
    Generation_, 1986, p 500. http://cg.scs.carleton.ca/~luc/rnbookindex.html
    """

    # p should be in (0.0, 1.0].
    if p <= 0.0 or p > 1.0:
    raise ValueError("p must be in the interval (0.0, 1.0]")
    elif p == 1.0:
    # If p is exactly 1.0, then the only possible generated value is 1.
    # Recognizing this case early means that we can avoid a log(0.0) later.
    # The exact floating point comparison should be fine. log(eps) works just
    # dandy.
    return 1

    # random() returns a number in [0, 1). The log() function does not
    # like 0.
    U = 1.0 - random.random()

    # Find the corresponding geometric variate by inverting the uniform variate.
    G = int(ceil(log(U) / log(1.0 - p)))
    return G

    --
    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, Jul 24, 2006
    #7
  8. MyInfoStation

    Paul Rubin Guest

    I usually owuld write that as int(ceil(log(U, 1.0 - p))).
     
    Paul Rubin, Jul 24, 2006
    #8
  9. MyInfoStation

    Robert Kern Guest

    Knock yourself out. I was cribbing from my C implementation in numpy.

    --
    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, Jul 24, 2006
    #9
  10. MyInfoStation

    Paul Rubin Guest

    Oh cool, I thought you were pasting from a Python implementation. No prob.
     
    Paul Rubin, Jul 24, 2006
    #10
    1. Advertisements

Ask a Question

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

You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.