基本的蒙特卡洛采样方法

# beta distribution with shape parameter a=3,b=6
a=3
b=6

# just show max value should be at 2.0/7 and max value is ceilled to 2 decimal points
beta_max = round(beta.pdf(2.0/7,a,b) + 0.005, 2)

# plot beta distribution
x = np.linspace(0.0,1.0,50000*beta_max)
plt.plot(x, beta.pdf(x, a, b),'r-', lw=1, alpha=0.6, label='beta pdf')
plt.plot(x,(np.ones(len(x)) * beta_max).tolist())
plt.show()

# random sampling 50000 on uniform in [0,1)
uni_samples = np.random.uniform(0.0,1.0,50000*beta_max)
# calculate "accept probability" = beta.pdf(x) / beta_max
# then do another uniform sample, if sample < beta.pdf(x) / beta_max,
# we accept the x as an sample from beta(3,6).
beta_samples = uni_samples[np.random.uniform(0.0,1.0,len(x)) < beta.pdf(uni_samples,a,b) / beta_max]

plt.hist(beta_samples, bins=50,normed=True)
plt.title("Reject Sampling for Beta(3,6) Histogram")
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def beta_pdf(x,a=3,b=6):
    return np.power(x,(a-1))*np.power((1-x),(b-1))

def f(x, a=3, b=6):
    """
    Log beta distribution
    There is constant has been remove which is gamma(a+b)/(gamma(a)+gamma(b))
    """
    return (a-1)*np.log(x)+(b-1)*np.log(1-x)

def fprima(x, a=3, b=6):
    """
    Derivative of Log beta distribution
    """
    return (a-1)/x-(b-1)/(1-x)

# 初始点
x = np.array([0.1,0.4,0.8])

# 初始化交点的横轴对应的值
z = np.zeros(len(x)+1)
z[0] = 0.0
z[-1] = 1.0

# 计算切线交点对应的横轴的值
# use np.diff function will give more simple code.
for j in range(len(x)-1):
  z[j+1] = (f(x[j+1])-f(x[j]) - x[j+1]*fprima(x[j+1]) + x[j]*fprima(x[j])) / (fprima(x[j]) - fprima(x[j+1]))

# 计算切线交点对应的纵轴的值
# h log pdf, hprime, the deriative of h
h = f(x)
hprime = fprima(x)

# the upper bound line function, which are piecewise liines,
# u(x) = h(x[j]) + (x-x[j])hprime(x[j])
# we only calculated the intersect points upper bound.
u = hprime[[0]+range(len(x))]*(z-x[[0]+range(len(x))]) + h[[0]+range(len(x))]

# 画出对应的pdf的对数曲线,切线,并标注交点
# plot log beta distribution
fig, ax = plt.subplots()
log_x = np.linspace(0.0,1.0,1000)
ax.plot(log_x,f(log_x))

# plot tangent lines
for i in range(len(x)):
    log_line_x = np.linspace(z[i],z[i+1],30)
    log_line_y = h[i] + hprime[i]*(log_line_x-x[i])
    ax.plot(log_line_x,log_line_y,color='green')

# mark x-axis points, i.e. x and z
xticks = []
xticklabels = []

for i in range(len(x)):    
    ax.plot([x[i]]*30,np.linspace(-35,0,30),ls='dotted',color='red')
    xticks += [x[i]]
    xticklabels += ["x[%d]"%i]

for i in range(len(z)):
    ax.plot([z[i]]*30,np.linspace(-35,0,30),ls='dotted',color='red')
    xticks += [z[i]]
    xticklabels += ["z[%d]"%i]

ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels)

# mark y-axis points, i.e. u
ax.scatter(z,u,color='blue')
fig.show()

s = np.hstack([0,np.cumsum(np.diff(np.exp(u))/hprime)])
cu = s[-1]
cdf=s/cu

def inverse_and_test(y):
    '''
    for an uniformed sample y, use inverse cdf get sample xt
    1. inverse piece-wised cdf to get piece
    2. on that piece, inverse cdf to get x
    '''     
    # Find the largest z such that cdf < y. i.e. inverse piecewise cdf get piece
    i = np.nonzero(cdf < y)[0][-1]

    # Figure out x from inverse cdf in relevant piece
    xt = x[i] + (-h[i] + np.log(hprime[i]*(cu*y - s[i]) + np.exp(u[i]))) / hprime[i]

    ut = h[i] + (xt-x[i])*hprime[i]
    ht = f(xt)
    if np.random.uniform() < np.exp(ht-ut):
       return xt

uni_samples = np.random.uniform(0.0,1.0,50000)
beta_samples = filter(lambda x: x is not None,map(inverse_and_test,uni_samples))

# plot sampling result
fig, ax = plt.subplots()
ax.hist(beta_samples, bins=50,normed=True)
fig.show()