Paul Rubin wrote: > Paul Rubin <http://[EMAIL PROTECTED]> writes: >> n = log(c, 1-p) - 1 > > I meant n = log(c/p, 1-p) - 1 > sorry.
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 -- http://mail.python.org/mailman/listinfo/python-list