马尔可夫链蒙特卡洛(Markov chain Monte Carlo/MCMC)学习

随机模拟的基本思想

假设有一个矩形区域S，其内部有一个不规则区域M，怎样求M的面积？

1. 将不规则区域M划分成很多规则的小区域，求所有小区域面积之和；
2. 随机的在S内部取点，计算点落在M内部的比例，该比例可近似S和M的面积比。

方法2采用的就是随机模拟的方法解决问题，随机模拟的关键在于如何将实际复杂的问题，转化为抽样可以解决的问题。

采样算法

给定统计样本集，如何估计产生这个样本集的随机变量概率密度函数，是我们比较熟悉的概率密度估计问题。求解概率密度估计问题的常用方法是最大似然估计、最大后验估计等。

但是，我们思考概率密度估计问题的逆问题:给定一个概率分布$p(x)$，如何让计算机生成满足这个概率分布的样本，这个问题就是统计模拟中研究的重要问题–采样(Sampling)。

计算机可以直接生成均匀分布的样本，现在已经有很多方法，那么我们的核心问题就是如何利用均匀抽样，借助随机模拟的方法，实现满足特定概率分布$p(x)$的抽样。

接受-拒绝抽样

对于非均匀分布，我们如果想要生成符合其分布的样本，就需要转换思路，利用均匀抽样和随机模拟来达到目的。

对于上图中的$\tilde{p}(z)$，是较为复杂的概率密度函数，我们通过以下步骤来获得符合其分布的样本：

1. 构造一个比较简单的容易抽样的高斯分布$kq(z)$，满足$kq(z)\ge \tilde{p}(z)$；
2. 对高斯分布进行采样，得到一个$z_0$；
3. 在$(z_0,kq(z_0))$之间进行均匀采样，得到${\mu }_0$；
4. 如果${\mu }_0>\tilde{p}(z_0)$，则拒绝该采样作为该分布的一个样本，反之则接受；
5. 重复2-4直到接受足够多的样本。

Python代码实现：

def rejection_sampling(iter=1000):
    samples = []

    for i in range(iter):
        z = np.random.normal(50, 30)
        u = np.random.uniform(0, k*q(z))

        if u <= p(z):
            samples.append(z)

    return np.array(samples)

这里有两个点需要理解，第一个点是$z_0$的采样，我们为什么要找一个和复杂分布类似的简单分布？目的就是在采样$z_0$的时候取到复杂分布取值范围内的概率比较大，避免过多的无用计算；第二个点是${\mu }_0$的采样，也是为了说明为什么拒绝采样是可行的，当$\tilde{p}(z_0)$较低时，决绝的概率较大，反之接受的概率较大，大致是符合$\tilde{p}(x)$的分布的。

这样，通过简单的高斯分布采样以及计算机可以实现的均匀采样，间接实现了复杂的分布的采样。

但是，这种方法在高维情况小的劣势也很明显，我们很难找到合适的$kq(x)$，不合适的上界造成的结果就是采样的拒绝率会超高，浪费了大量的计算。

马尔科夫链及其平稳分布

马尔科夫链由一个状态空间下的初始状态概率 π = { π 0 , π 1 , . . . , π n } \pi=\{\pi_0, \pi_1, ..., \pi_n\} π={ π0​,π1​,...,πn​},和状态转移矩阵$P$组成，${\pi }_i$表示初始状态是状态$i$的概率，$P_{ij}$表示状态$i$转换为状态$j$的概率，根据初始状态概率和状态转移矩阵，我们可以得到一个状态序列 X = { X 1 , X 2 , . . . , X n } X=\{X_1, X_2, ..., X_n\} X={ X1​,X2​,...,Xn​}：
$$P(X_{t+1}=x|X_t,X_{t-}$$