In [31]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sps
import math
"""
kg = -22.75
Ng = +21.26
Ag = 1.27
pdf = kg* (N/Ng)**Ag*math.exp(-N/Ng)
-22.75*(x/21.26)^(-1.27)*exp(-x/21.26)
gamma_pdf(x)=x^(k-1)*exp(x/b)/(b^k*gamma(k))
"""
shape, scale =3, 2.
s = np.random.gamma(shape, scale, 1000)
count, bins, ignored = plt.hist(s, 50, normed=True)
y = bins**(shape-1)*(np.exp(-bins/scale) / (sps.gamma(shape)*scale**shape))
plt.plot(bins, y, linewidth=2, color='r')
plt.show()
In [33]:
t = np.arange(15., 24., 0.2)
#plt.plot(t, t, 'r--', t, t**2, 'bs', t, t**3, 'g^')
plt.plot(t,-22.75*(t/21.26)**(-1.27)*np.exp(-t/21.26))
plt.show()
In [24]:
class GeneralRandom:
"""This class enables us to generate random numbers with an arbitrary
distribution."""
def __init__(self, x = np.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 = np.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 = np.arange(Nrl)/float(Nrl)
delta = 1.0/Nrl
self.inversecdf = np.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 = np.concatenate((np.diff(self.inversecdf), [0]))
def random(self, N = 1000):
"""Give us N random numbers with the requested distribution"""
idx_f = np.random.uniform(size = N, high = self.Nrl-1)
idx = np.array([idx_f],'i')
y = self.inversecdf[idx] + (idx_f - idx)*self.delta_inversecdf[idx]
return y
def plot_pdf(self):
plt.plot(self.x, self.pdf)
def self_test(self, N = 1000):
plt.figure()
#The cdf
plt.subplot(2,2,1)
plt.plot(self.x, self.cdf)
#The inverse cdf
plt.subplot(2,2,2)
y = np.arange(self.Nrl)/float(self.Nrl)
plt.plot(y, self.inversecdf)
#The actual generated numbers
plt.subplot(2,2,3)
y = self.random(N)
p1, edges = np.histogram(y, bins = 50,
range = (self.x.min(), self.x.max()),
normed = True, new = True)
x1 = 0.5*(edges[0:-1] + edges[1:])
plt.plot(x1, p1/p1.max())
plt.plot(self.x, self.pdf/self.pdf.max())
In [30]:
# Generating N random numbers that probability distribution
# fits to any given function curve
# FB - 201006137
import math
import random
# define any function here!
def f(x):
return -22.75*(x/21.26)**(-1.27)*math.exp(-x/21.26)
# f(x) = 1.0 : for uniform probability distribution
# f(x) = x : for triangular probability distribution
# (math.sqrt(random.random()) would also produce triangular p.d. though.)
# f(x) = math.exp(-x*x/2.0)/math.sqrt(2.0*math.pi) : for std normal p.d.
# (taking average of (last) 2,3,... random.random() values would also
# produce normal probability distributions though.)
# define any xmin-xmax interval here! (xmin < xmax)
xmin = 14.0
xmax = 25.0
# find ymin-ymax
numSteps = 1000000 # bigger the better but slower!
ymin = f(xmin)
ymax = ymin
for i in range(numSteps):
x = xmin + (xmax - xmin) * float(i) / numSteps
y = f(x)
if y < ymin: ymin = y
if y > ymax: ymax = y
n = 10 # how many random numbers to generate
for i in range(n):
while True:
# generate a random number between 0 to 1
xr = random.random()
yr = random.random()
x = xmin + (xmax - xmin) * xr
y = ymin + (ymax - ymin) * yr
if y <= f(x):
print xr
break
"""
s = np.random.gamma(shape, scale, 1000)
count, bins, ignored = plt.hist(s, 50, normed=True)
y = bins**(shape-1)*(np.exp(-bins/scale) / (sps.gamma(shape)*scale**shape))
plt.plot(bins, y, linewidth=2, color='r')
plt.show()
"""
Out[30]:
In [ ]: