Welcome, guest | Sign In | My Account | Store | Cart

This is a class that allows you to set up an arbitrary probability distribution function and generate random numbers that follow that arbitrary distribution.

Python, 75 lines
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75``` ```import pylab import numpy class GeneralRandom: """This class enables us to generate random numbers with an arbitrary distribution.""" def __init__(self, x = pylab.arange(-1.0, 1.0, .01), p = None, Nrl = 1000): """Initialize the lookup table (with default values if necessary) Inputs: x = random number values p = probability density profile at that point Nrl = number of reverse look up values between 0 and 1""" if p == None: p = pylab.exp(-10*x**2.0) self.set_pdf(x, p, Nrl) def set_pdf(self, x, p, Nrl = 1000): """Generate the lookup tables. x is the value of the random variate pdf is its probability density cdf is the cumulative pdf inversecdf is the inverse look up table """ self.x = x self.pdf = p/p.sum() #normalize it self.cdf = self.pdf.cumsum() self.inversecdfbins = Nrl self.Nrl = Nrl y = pylab.arange(Nrl)/float(Nrl) delta = 1.0/Nrl self.inversecdf = pylab.zeros(Nrl) self.inversecdf[0] = self.x[0] cdf_idx = 0 for n in xrange(1,self.inversecdfbins): while self.cdf[cdf_idx] < y[n] and cdf_idx < Nrl: cdf_idx += 1 self.inversecdf[n] = self.x[cdf_idx-1] + (self.x[cdf_idx] - self.x[cdf_idx-1]) * (y[n] - self.cdf[cdf_idx-1])/(self.cdf[cdf_idx] - self.cdf[cdf_idx-1]) if cdf_idx >= Nrl: break self.delta_inversecdf = pylab.concatenate((pylab.diff(self.inversecdf), [0])) def random(self, N = 1000): """Give us N random numbers with the requested distribution""" idx_f = numpy.random.uniform(size = N, high = self.Nrl-1) idx = pylab.array([idx_f],'i') y = self.inversecdf[idx] + (idx_f - idx)*self.delta_inversecdf[idx] return y def plot_pdf(self): pylab.plot(self.x, self.pdf) def self_test(self, N = 1000): pylab.figure() #The cdf pylab.subplot(2,2,1) pylab.plot(self.x, self.cdf) #The inverse cdf pylab.subplot(2,2,2) y = pylab.arange(self.Nrl)/float(self.Nrl) pylab.plot(y, self.inversecdf) #The actual generated numbers pylab.subplot(2,2,3) y = self.random(N) p1, edges = pylab.histogram(y, bins = 50, range = (self.x.min(), self.x.max()), normed = True, new = True) x1 = 0.5*(edges[0:-1] + edges[1:]) pylab.plot(x1, p1/p1.max()) pylab.plot(self.x, self.pdf/self.pdf.max()) ```

#### 1 comment

joel.thorarinson 14 years, 10 months ago

Hey awesome little program, why did someone downvote it?

 Created by Kaushik Ghose on Wed, 5 Nov 2008 (MIT)