一. MCMC算法定义

    MCMC是Markov chain Monte Carlo的缩写，即马尔可夫链蒙特卡罗方法。MCMC是一组利用马尔可夫链从随机分布中取样的算法。生成的马氏链即是对于目标分布的近似估计。

常见算法：

Metropolis-Hastings算法
Gibbs取样算法
Hamiltonian Monte Carlo算法（效率最高，但不常用）
在统计推断和贝叶斯推断中，常用的是前两种算法。

二、Metropolis算法

假设随机变量X服从一个概率密度函数P。

初始化X，即给定随机变量一个起始点 $X_{old}$
生成随机的跳跃 $\Delta X$ ~ $N(0,\sigma)$ ， $X_{new} = X_{old}+\Delta X$
计算随机变量移动到提议的点的概率
$P_{moving} = min(1, \frac{X_{new}} {X_{old}})$
生成一个随机数 $C$ ~ $U (0, 1)$ ，如果 $C < P_{moving}$ 则接受跳跃，随机变量取新的值 $X_{new}$ ；反之则任取原来的值。
重复步骤2~4，并记录下随机变量所有的取值，知道蒙特卡洛模拟生成的随机游走足以代表目标分布。

利用Gamma分布对Metropolis算法进行测试

# Metropolis算法
import numpy
from scipy.special import gamma
import scipy.stats
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def MCMC(P,X0,chain,space):
    '''
    P - 随机变量X服从的概率密度函数
    X0 - MCMC的起始值
    chain - MCMC的长度
    space - 随机变量的取值范围，如[0,inf]
    '''
    
    if not callable(P):
        raise Exception('P必须是函数')
        
    X_current = X0
    X = [X_current]
    
    while True:
        Delta_X = scipy.stats.norm(loc=0,scale=2).rvs()
        X_proposed = X_current + Delta_X
        
        if X_proposed < space[0] or X_proposed > space[1]:
            p_moving = 0
        elif P(X_current) == 0:
            p_moving = 1
        else:
            p_moving = min(1,P(X_proposed)/P(X_current))
            
        if scipy.stats.uniform().rvs() <= p_moving:  #新生成的服从0，1均匀分布的值小于等于p_moving
            X.append(X_proposed)
            X_current = X_proposed
        else:
            X.append(X_current)
            
        if len(X) >= chain:
            break
        
    return numpy.array(X)

# 测试Metropolis算法：生成服从Gamma分布的马氏链
def GammaDist(x,k,theta):
    '''定义Gamma分布'''
    return 1/(gamma(k)*theta**k) * x**(k-1) * numpy.exp(-x/theta)

# 得到参数k=2，theta=2的分布
gammadist = lambda x: GammaDist(x,2,2)

X = MCMC(gammadist,2,100,[0,numpy.inf])   #最终结果：有MCMC算法生成的马氏链

这里生成的X就是使用MCMC算法生成的近似服从Gamma分布，是对目标分布的近似估计。下面通过动画演示这一过程。

# 用一个动画来演示生成的马氏链
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))
fig.suptitle('MCMC Process')


def anima_chain(index):
    ax1.clear()

    ax1.plot(X[:index + 1], 'r-')
    # ax1.set_xlim([1,1000])
    ax1.set_xlabel('chain')  # 第多少步
    ax2.set_label('X')  # 生成的马氏链

    # ax1.set_title('MCMC Process')


def anima_density(index):
    ax2.clear()

    x = numpy.arange(0., 12, .1)
    ax2.plot(x, gammadist(x), 'k-')
    if index <= 10:
        y = numpy.zeros(len(x))
        ax2.plot(x, y)
    else:
        density = scipy.stats.gaussian_kde(X[:index - 1])
        ax2.plot(x, density(x), 'r-')

    ax2.set_xlim([0.,12])
    ax2.set_ylim([0,.5])
    ax2.set_xlabel('X')
    ax2.set_ylabel('density')

    # ax2.set_title('MCMC density estimate')

ani_chain = animation.FuncAnimation(fig, anima_chain, interval=1)
ani_density = animation.FuncAnimation(fig, anima_density, interval=1)


def main():
    plt.show()

if __name__ == '__main__':
    main()