马尔科夫蒙特卡洛_吉布斯抽样算法(Markov Chain Monte Carlo(MCMC)_Gibbs Sampling)

定义

输入:目标概率分布的密度函数$p(x)$,函数$f(x)$

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

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

(1)初始化。给出初始样本$x^0=(x_1^{(0)},x_2^{(0)},\cdots ,x_k^{(0)}{)}^T$。

(2)对i循环执行

设第$(i-1)$次迭代结束时的样本为$x^{(i-1)}=(x_1^{(i-1)},x_2^{(i-1)},\cdots ,x_k^{(i-1)}{)}^T$,则第$i$次迭代进行如下几步操作:

$\{\begin{array}{l}(1)由满足条件分布p(x_1|x_2^{(i)},\cdots ,x_k^{(i-1)})抽取x_1^{(i)}\\ ⋮\\ (j)由满足条件分布p(x_j|x_1^{(i)},\cdots ,x_{j-1}^{(i)},x_{j+1}^{(i)},\cdots ,x_k^{(i-1)})抽取x_j^{(i)}\\ ⋮\\ (k)由满足条件分布p(x_k|x_1^{(i)},\cdots ,x_{(k-1)}^{(}i)抽取x_k^{(i)}\end{array}$

得到第$i$次迭代值$x^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots ,x_k^{(i)}{)}^T$
(3)得到样本集合
{ x ( m + 1 ) , x ( m + 2 ) , ⋯   , x ( n ) } \{ x^{(m+1)},x^{(m+2)},\cdots,x^{(n)} \} {x(m+1),x(m+2),⋯,x(n)}
(4)计算
$$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 gibbs_sampling(dim, conditional_sampler, x0=None, burning_steps=1000, max_steps=10000, epsilon=1e-8, verbose=False):
    """
    给定一个从p（x_j|x_1，x_2，…x_n）采样的条件采样器
    返回样本x~p的列表，其中p是条件分布的原始分布。
    x0是x的初始值。如果未指定，则将其设置为零向量。
    条件采样器将（x，j）作为参数
    """
    x = np.zeros(dim) if x0 is None else x0
    samples = np.zeros([max_steps - burning_steps, dim])
    for i in range(max_steps):
        for j in range(dim):
            x[j]  = conditional_sampler(x, j)
            if verbose:
                print("New value of x is", x_new)
        if i >= burning_steps:
            samples[i - burning_steps] = x
    return samples


if __name__ == '__main__':
    i = 0
    def demonstrate(dim, p, desc, **args):
        samples = gibbs_sampling(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.savefig(str(i)+".png")
        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 gaussian_sampler1(x, j):
        return np.random.normal()
    demonstrate(2, gaussian_sampler1, "Gaussian distribution with mean of 2 and 3")

    # example 2:
    mean = np.array([2, 3])
    covariance = np.array([[1, 0],
                           [0, 1]])
    covariance_inv = np.linalg.inv(covariance)
    det_convariance = 1
    def gaussian_sampler2(x, j):
        if j == 0:
            return np.random.normal(2)
        else:
            return np.random.normal(3)
    demonstrate(2, gaussian_sampler2, "Gaussian distribution with mean of 2 and 3")

    # example 3:
    def blocks_sampler(x, j):
        sample = np.random.random()
        if sample > .5:
            sample += 1.
        return sample
    demonstrate(2, blocks_sampler, "Four blocks")

    # example 4:
    def blocks_sampler(x, j):
        sample = np.random.random()
        if sample > .5:
            sample += 100.
        return sample
    demonstrate(2, blocks_sampler, "Four blocks with large gap.")