马尔科夫蒙特卡洛_梅特罗波利斯－黑斯廷斯算法(Markov Chain Monte Carlo(MCMC)_Metropolis Hastings algorithm)

定义

输入:抽样的目标分布的密度函数$p(x)$,函数$f(x)$;

输出:$p(x)$的随机样本$x_{m+1},x_{m+2},\cdots ,x_n$,函数样本均值$f_{mn};$

参数:收敛步数$m$,迭代步数$n$。

(1)任意选择一个初始值$x_0$

(2)对$i=1,2,\cdots ,n$循环执行

      (a)设状态$x_{i-1}=x,按照建议分布g(x,x^{{}^{\prime }})随机抽取一个候选状态x^{{}^{\prime }}$。

      (b)计算接受概率
α ( x , x ′ ) = m i n { 1 , p ( x ′ ) q ( x ′ , x ) p ( x ) q ( x , x ′ ) } \alpha(x,x^{'}) = min \{ 1,\frac{p(x^{'})q(x^{'},x)}{p(x)q(x,x^{'})} \} α(x,x′)=min{1,p(x)q(x,x′)p(x′)q(x′,x)​}
      ©从区间$(0,1)$中按均匀分布随机抽取一个数$\mu$。若$\mu \le \alpha (x,x^{{}^{\prime }}),则状态x_i=x^{{}^{\prime }};否则,状态x_i=x$。

(3)得到样本集合 { x m + 1 , x m + 2 , ⋯   , x n } \{x_{m+1},x_{m+2},\cdots,x_n\} {xm+1​,xm+2​,⋯,xn​}

计算
$$f_{mn}=\frac{1}{n-m}\sum _{i=m+1}^{n}f(x_i)$$

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import gaussian_kde


def gaussian_kernel(x1, x2):
    return np.exp(-((x1 - x2) ** 2).sum())

def gaussian_sampler(x):
    return np.random.normal(x)

def metropolis_hasting(dim, p, q=gaussian_kernel, q_sampler=gaussian_sampler, x0=None, burning_steps=1000, max_steps=10000, epsilon=1e-8, verbose=False):
    """
    给定分布函数p（它不需要是概率，似然函数就足够了），
    以及推荐分布q，
    返回样本x~p的列表，
    其中样本数为max_steps-burning_steps。
    q_sampler是一个以x为输入并返回q（x_new|x_old）样本的函数。
    q是表示q（x_new|x_old）的分布函数。
    q取（xold，xnew）作为参数。
    """
    x = np.zeros(dim) if x0 is None else x0
    samples = np.zeros([max_steps - burning_steps, dim])
    for i in range(max_steps):
        x_new = q_sampler(x)
        accept_prob = (p(x_new) + epsilon) / (p(x) + epsilon) * q(x, x_new) / q(x_new, x)
        if verbose:
            print("New value of x is", x_new)
        if np.random.random() < accept_prob:
            x = x_new
        elif verbose:
            print("New value is dropped")
        if i >= burning_steps:
            samples[i - burning_steps] = x
    return samples


if __name__ == '__main__':
    def demonstrate(dim, p, desc, **args):
        samples = metropolis_hasting(dim, p, **args)
        z = gaussian_kde(samples.T)(samples.T)
        plt.scatter(samples[:, 0], samples[:, 1], c=z, marker='.')
        plt.plot(samples[: 100, 0], samples[: 100, 1], 'r-')
        plt.title(desc)
        plt.show()

    # example 1:
    mean = np.array([2, 3])
    covariance = np.array([[1, 0],
                           [0, 1]])
    covariance_inv = np.linalg.inv(covariance)
    det_convariance = 1
    def gaussian1(x):
        return np.exp(-.5 * (x - mean).T @ covariance_inv @ (x - mean))
    demonstrate(2, gaussian1, "Gaussian distribution with mean of 2 and 3")

    # example 2:
    mean = np.array([2, 3])
    covariance = np.array([[1, .5],
                           [.5, 1]])
    covariance_inv = np.linalg.inv(covariance)
    det_convariance = 1
    def gaussian2(x):
        return np.exp(-.5 * (x - mean).T @ covariance_inv @ (x - mean))
    demonstrate(2, gaussian2, "Gaussian distribution with mean of 2 and 3")

    # example 3:
    def blocks(x):
        if (0 < x[0] < 1 or 2 < x[0] < 3) and (0 < x[1] < 1 or 2 < x[1] < 3):
            return 1
        return 0
    demonstrate(2, blocks, "Four blocks")

    # example 4:
    def blocks(x):
        if (0 < x[0] < 1 or 200 < x[0] < 300) and (0 < x[1] < 1 or 200 < x[1] < 300):
            return 1
        return 0
    demonstrate(2, blocks, "Four blocks with large gap. (Monte Carlo doesn't solve everything)")