import numpy as np
import matplotlib.pyplot as plt
# probability distribution we're trying to calculate
p = lambda x: np.exp(-x)
# CDF of p
CDF = lambda x: 1-np.exp(-x)
# invert the CDF
invCDF = lambda x: -np.log(1-x)
# domain limits
xmin = 0 # the lower limit of our domain
xmax = 6 # the upper limit of our domain
# range limits
rmin = CDF(xmin)
rmax = CDF(xmax)
N = 10000 # the total of samples we wish to generate
# generate uniform samples in our range then invert the CDF
# to get samples of our target distribution
R = np.random.uniform(rmin, rmax, N)
X = invCDF(R)
# get the histogram info
hinfo = np.histogram(X,100)
# plot the histogram
plt.hist(X,bins=100, label=u'Samples');
# plot our (normalized) function
xvals=np.linspace(xmin, xmax, 1000)
plt.plot(xvals, hinfo[0][0]*p(xvals), 'r', label=u'p(x)')
# turn on the legend
plt.legend()
plt.show()

# -*- coding: utf-8 -*-
'''
The following code produces samples that follow the distribution P(x)=e^−x 
for x=[0, 10] and generates a histogram of the sampled distribution.
'''
import numpy as np
import matplotlib.pyplot as plt
P = lambda x: np.exp(-x)
# domain limits
xmin = 0 # the lower limit of our domain
xmax = 10 # the upper limit of our domain
# range limit (supremum) for y
ymax = 1
N = 10000  # the total of samples we wish to generate
accepted = 0 # the number of accepted samples
samples = np.zeros(N)
count = 0  # the total count of proposals
# generation loop
while (accepted < N):
    # pick a uniform number on [xmin, xmax) (e.g. 0...10)
  x = np.random.uniform(xmin, xmax)
    # pick a uniform number on [0, ymax)
  y = np.random.uniform(0,ymax)
    # Do the accept/reject comparison
  if y < P(x):
    samples[accepted] = x
    accepted += 1
    count +=1
  print count, accepted
# get the histogram info
# If bins is an int, it defines the number of equal-width bins in the given range 
(n, bins)= np.histogram(samples, bins=30) # Returns: n-The values of the histogram,n是直方圖中柱子的高度
# plot the histogram
plt.hist(samples,bins=30,label=u'Samples')  # bins=30即直方圖中有30根柱子
# plot our (normalized) function
xvals=np.linspace(xmin, xmax, 1000)
plt.plot(xvals, n[0]*P(xvals), 'r', label=u'P(x)')
# turn on the legend
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
p = lambda x: np.exp(-x)     # our distribution
g = lambda x: 1/(x+1)      # our proposal pdf (we're choosing k to be 1)
CDFg = lambda x: np.log(x +1)  # generates our proposal using inverse sampling
# domain limits
xmin = 0 # the lower limit of our domain
xmax = 10 # the upper limit of our domain
# range limits for inverse sampling
umin = CDFg(xmin)
umax = CDFg(xmax)
N = 10000 # the total of samples we wish to generate
accepted = 0 # the number of accepted samples
samples = np.zeros(N)
count = 0 # the total count of proposals
# generation loop
while (accepted < N):
    # Sample from g using inverse sampling
  u = np.random.uniform(umin, umax)
  xproposal = np.exp(u) - 1
  # pick a uniform number on [0, 1)
  y = np.random.uniform(0, 1)
  # Do the accept/reject comparison
  if y < p(xproposal)/g(xproposal):
    samples[accepted] = xproposal
    accepted += 1
    count +=1
  print count, accepted
# get the histogram info
hinfo = np.histogram(samples,50)
# plot the histogram
plt.hist(samples,bins=50, label=u'Samples');
# plot our (normalized) function
xvals=np.linspace(xmin, xmax, 1000)
plt.plot(xvals, hinfo[0][0]*p(xvals), 'r', label=u'p(x)')
# turn on the legend
plt.legend()
plt.show()