Inverse Transform Sampling

ITS.py

"""
Author: Sen MU
Inverse transform sampling
reference article: http://usmanwardag.github.io/python/astronomy/2016/07/10/inverse-transform-sampling-with-python.html
"""
"""
steps:
1- Normalize a distribution in terms of its CDF (cumulative distribution function).

2- Generate a random number u from standard uniform distribution in interval [0, 1].

3- Compute an event x from the distrubtion such that f(x) = u.

4- Take x to be the random event drawn from the distribtion.
"""

import numpy as np
from matplotlib import pyplot as plt
import numpy.random as random

# Step 1: Generate Cumulative Distribution Function
data = [[1, 2, 3, 4, 5, 6, 7, 8, 9],[2000, 4040, 6500,1000,4000, 6000, 5000, 4020, 2070]]
x = np.array(data[0])
y = np.array(data[1])

# calculate cumulative probability
prob = y/float(sum(y))
cumulative_prob = np.cumsum(prob)

# Step 2: Generate Random Numbers
N = 5000  # N is the amount of random numbers
random_numbers = random.uniform(0,1,N)

# Step 3: Generate Samples using inverse function of CDF
gen_x = []
for r in random_numbers:
    gen_x.append(int(x[np.argwhere(cumulative_prob == min(cumulative_prob[(cumulative_prob - r) > 0]))]))

# Step 4: Test the sampling result
gen_x = ((np.array(gen_x) - 1) / 1).astype(int)
times = np.arange(1, y.shape[0]+1, 1) 
lc = np.bincount(gen_x, minlength=len(times))

print('generated: ', lc/float(sum(lc)))
print('origin: ', prob)

plot1, = plt.plot(lc/float(sum(lc)), 'r', label='Generated distribution')
plot2, = plt.plot(prob,'g',label='Original distribution')
plt.xlabel('x')
plt.ylabel('Probability')
plt.legend(handles=[plot1,plot2])
plt.show()