MCMC（Markov Chain Monte Carlo）的理解与实践（Python）

Markov Chain Monte Carlo (MCMC) methods are a class of 
algorithms for sampling from a probability distribution
based on constructing a Markov chain that has the desired 
distribution as its stationary distribution. The state of 
the chain after a number of steps is then used as a sample 
of the desired distribution. The quality of the sample 
improves as a function of the number of steps.


With MCMC, we draw samples from a (simple) proposal 
distribution so that each draw depends only on the state of 
the previous draw (i.e. the samples form a Markov chain,
 θp=θ+Δθ,Δθ∼N(0,σ)θp=θ+Δθ,Δθ∼N(0,σ)).

import numpy as np
import scipy.special as ss
import matplotlib.pyplot as plt

def beta_s(x, a, b):
    return x**(a-1)*(1-x)**(b-1)
def beta(x, a, b):
    return beta_s(x, a, b)/ss.beta(a, b)

def plot_mcmc(a, b):
    cur = np.random.rand()
    states = [cur]
    for i in range(10**5):
        next, u = np.random.rand()
        if u < np.min((beta_s(next, a, b)/beta_s(cur, a, b), 1)):
            states.append(next)
            cur = next
    x = np.arange(0, 1, .01)
    plt.figure(figsize=(10, 5))
    plt.plot(x, beta(x, a, b), lw=2, label='real dist: a={}, b={}'.format(a, b))
    plt.hist(states[-1000:], 25, normed=True, label='simu mcmc: a={}, b={}'.format(a, b))
    plt.show()

if __name__ == '__main__':
    plot_mcmc(0.1, 0.1)
    plot_mcmc(1, 1)
    plot_mcmc(2, 3)