MCMC抽样法（学习笔记）

前言：

后验分布种给出了每个参数值的可能性，但很难看出完整的分布，也无法用分析解决。于是我们需要MCMC法；
MCMC法有两大要素：蒙特卡罗模型和马尔可夫链

蒙特卡罗模型

蒙特卡罗模型：重复生成（大量）随机数来估计固定参数（仅给出一个参数的近似值但不能直接计算）
- 举例：估算不规则区域面积
  
  image.png
通过随机生成20个点来估计圆面积（ $\approx \frac{15}{20}*100=75$ ）
应用：
1. 估算不规则区域面积
2. 模拟复杂过程（预测天气、估计选举结果......)

马尔可夫链

历史发展：

image.png
- 如图，钟形曲线被看作自然界中常见模式，但钟形线甚至大多数规律所描述的一系列事件，彼此都是完全独立的。因此，俄罗斯数学家和神学家帕维尔·涅克拉索夫（Pavel Nekrasov）认为：现实世界中的事物相互依存，不符合数学模型或分布！
- 但是安德烈·马尔科夫（Andrey Markov）试图证明非独立事件也可能符合数学模型。马尔可夫链一个著名的应用案例是——从一部俄罗斯诗歌作品中输出数千个“两字母对”（two-character pairs）。利用这些对，计算每一字母的条件概率。也就是说，给定前面的字母或空格，就能预测下一个字母是A还是T，或是一个空格。利用这些概率，马尔可夫可以模拟任意长的字符序列，这就是马尔可夫链。当字符串一定长的时候，起始字母的选择对后续字母的影响可以忽略不计即字母的分布遵循一种模式。也就是说，相互依赖的时间，如果它们受到固定概率的影响，也将遵循平均水平模式！
马尔可夫链有概率相关的序列组成，每一事件来自一组结果，每个结果根据一组固定的概率确定下一个结果。
一种重要特征——无记忆：预测下一个时间所需要的所有信息都能从当前状态找到，从历史事件中找不到新的信息（类“蛇梯棋”）
- 数学定义：
$p(E_i|E_i,...,E_1)=p(E_i|E_{i-1})$
- 一个形象的比喻：假如每天的天气是一个状态的话，那今天是不是晴天只依赖于昨天的天气而和前天的天气没有任何关系！
- 非周期的马尔可夫链具有收敛性质：如果一个非周期的马尔可夫链有状态转移矩阵P，且它的任何两个状态是连通的，那么转移n次后的状态与初始状态无关！
  - 例^[1]：对于如下图的股市模型：
    
    image.png
    
    共有三种状态，定义状态概率转化矩阵为P，且P(i,j)=P(j|i)(从状态i转化到状态j的概率)：
    
    $P= \begin{pmatrix} 0.9&0.075&0.025\\ 0.15&0.8&0.05\\ 0.25&0.25&0.5 \end{pmatrix}$
    我们会发现，对于任意初始的股市概率，假设初始概率分布：[0.3,0.4,0.3],代入转移矩阵计算状态2，3，……，n次状态转移（本例中约为6次）会发现股市概率分布稳定在：[0.625,0.3125,0.0625]这组概率分布上（不论初始股市概率是什么）。
- 利用MCMC模型估计后验分布形状
  
  利用MCMC方法,有效从后验分布种抽取样本，然后计算样本中的平均值
  - 步骤：
    
    输入选定的马尔可夫链状态转移矩阵Q，设定状态转移次数阈值n₁ ，需要的样本个数n₂ $\longrightarrow$ 从任意简单概率分布取一个初始状态值x₀（或自己给定） $\longrightarrow$ 重复X步骤n₁+n₂ -1次 $\longrightarrow$ 获得需要的平稳分布对应的样本集( $x_{n1}, x_{n_1+1},...,x_{n_!+n_2+-1}$ ) $\longrightarrow$ 得到概率分布
    - X步骤：(t表示重复的次数):
      1. 从条件概率分布Q(x|x_t)中取得样本x_*
      2. 均匀分布采样u~uniform[0,1] （获得一个随机数u）
      3. 如果u<α（x_t ,x_*）(接受率，α（x_t ,x~~）= $\pi$ (x~_)Q(x~, x_t)；若为M-H分布，则α(i,j)=min{ $\frac{\pi (j)Q(j,i)}{\pi (i)Q(i,j)}$ ,1},提高了采样效率)，则接受转移x_t $\rightarrow x_*$ ，即x_t+1 =x~~
      4. 否则不接受转移，即 $x_{t+1}=x_t$

利用MCMC模型估计后验分布形状

利用MCMC方法,有效从后验分布种抽取样本，然后计算样本中的平均值

步骤：

输入选定的马尔可夫链状态转移矩阵Q，设定状态转移次数阈值n₁ ，需要的样本个数n₂ $\longrightarrow$ 从任意简单概率分布取一个初始状态值x₀（或自己给定） $\longrightarrow$ 重复X步骤n₁+n₂ -1次
- X步骤：(t表示重复的次数)
1. 从条件概率分布Q(x|x_t)中取得样本x_*
2. 均匀分布采样u~uniform[0,1] （获得一个随机数u）
3. 如果u<α（x_t ,x_*）(接受率，α（x_t ,x~~）= $\pi$ (x~_)Q(x~, x_t))，则接受转移x_t $\rightarrow x_*$ ，即x_t+1 =x~~
4. 否则不接受转移，即 $x_{t+1}=x_t$
$\longrightarrow$ 获得需要的平稳分布对应的样本集( $x_{n1}, x_{n_1+1},...,x_{n_!+n_2+-1}$ ) $\longrightarrow$ 得到概率分布
例：MCMC单参数运用——估计人类平均身高
1. 人类平均身高的后验分布同时受到观测曲线和先验分布的影响，假设人类平均身高的后验分布如图：
  
  image.png
2. 对单参数（x轴上的数值），MCMC法沿x轴开始随机采样，在这里后验概率影响Q，则：
  
  image.png
  
  （集中的点垂直叠加，便于观察，事实上应该都是取在x轴上的）
3. 做直方图，获取概率分布：
  
  image.png

REFERENCE

《告别数学公式，图文解读什么是马尔可夫链蒙特卡罗方法（MCMC）》：https://zhuanlan.zhihu.com/p/32982140
《MCMC(三)MCMC采样和M-H采样》：https://www.cnblogs.com/pinard/p/6638955.html

来自wiki ↩