In [9]:

    

import matplotlib.pyplot as plt
import math
import random

# define any function here!
def f(x):
    log10kg =-22.75
    log10Ng = 21.26
    Ag = -1.27
    return log10kg + Ag*(x-log10Ng) - 10**x/10**log10Ng*math.log10(math.exp(1))

# 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 = 20.3
xmax = 22.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


xdata = []        
n = 1000 # 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):
            xdata.append(x)
            break

plt.hist(xdata,bins=25)
plt.title("Custom Histogram")
plt.xlabel("Value")
plt.ylabel("Frequency")
plt.show()            
            
"""
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[9]:





"\ns = np.random.gamma(shape, scale, 1000)\ncount, bins, ignored = plt.hist(s, 50, normed=True)\ny = bins**(shape-1)*(np.exp(-bins/scale) / (sps.gamma(shape)*scale**shape))\nplt.plot(bins, y, linewidth=2, color='r')\nplt.show()\n"

In [ ]: