5.1 马尔可夫链蒙特卡洛采样

本节预览：

5.1.1 蒙特卡洛

什么是蒙特卡洛？

蒙特卡洛（Monte Carlo）是一种思想或一类方法，它通过随机采样来估计数学问题的数值解，本质上是通过大量采样来暴力解决问题，特别适用于复杂的积分、概率计算和优化问题等。

例子：如何利用蒙特卡洛方法估计圆周率 $\pi$ ？

问题描述：在一个边长为 2 的正方形中，内嵌一个半径为 1 的圆，我们希望估计圆的面积。
生成随机样本：在正方形内随机生成 $n$ 个点，记录落在圆内的点数 $m$
计算估计值：
- 圆的面积与正方形面积的比值是 $\frac{\pi r^2}{(2r)^2} = \frac{\pi}{4}$
- 通过 $\frac{m}{n}$ 估计这个比值， $\pi \approx 4 \times \frac{m}{n}$
增加样本数量：随着样本数量 $n$ 增加， $\pi$ 的估计值会越来越接近真实值

在未知 $\pi$ 的情况下，我们无法直接计算圆的面积，于是通过暴力采样和计算点数量的方式估计圆的面积，进而估计 $\pi$ 的值。

蒙特卡洛如何帮助我们处理复杂后验分布？

在求解后验分布的时候，我们经常要在分布上做一些计算，比如求期望、方差、或其他积分。但是如果碰上复杂分布、高维变量，就会让计算变得特别困难。这时候我们就可以尝试绕过解析解，通过大量采样的方法直接得到统计量。

例如，我们想对一个连续随机变量 $X$ 求期望，概率分布函数为：

p(x) = e^{-0.5 \cdot x^2} + \sin(x) + x^3

这是一个非标准概率密度函数，我们无法直接求积分，但是我们可以从 $p(x)$ 中采样 $x_1, x_2, x_3, ...x_i$ ，再利用大数定理（在本书中不做详细解释）近似期望：

\mathbb{E}[X] = \int x_i\cdot p(x_i)\,d(x)\approx \frac{1}{n} \sum_{i}^{n} x_i

同样，对离散连续变量 $X$ 求期望，虽然不涉及积分困难，但当变量的维度很高的时候 $(X \in \mathbb{R}^D )$ ，会出现维数灾难，求和计算难度随着维度增加而指数级增加。我们也可以通过从 $P(x)$ 采样来近似期望：

\mathbb{E}[X] = \sum_{i}^{n} x_i\cdot P(x_i)\approx \frac{1}{n} \sum_{i}^{n}x_i

一句话总结：遇到棘手的概率分布形式，避开难求的解析解，通过大量采样进行暴力计算估计。

注意：这里我们已知概率函数 $p(x)$ , 只是这个 $p(x)$ 可能比较复杂，直接积分比较难求。

蒙特卡洛的实际应用案例

接下来我们借着一个模拟股票走势的代码案例，展示蒙特卡洛的用法。

假定某个股票目前的股价是 $S_t$ , 预测时间 $\Delta t$ 之后的股价 $S_{t+1}$ ，金融工程中有公式：

S_{t+1} = S_t + \hat{\mu} S_t \Delta t + \sigma S_t \epsilon \sqrt{\Delta t}

假定初始股价 $S_0=10$ ，我们可以模拟一下，一年之后 (244 个交易日 ) 的股价分布。

首先，定义初始状态和模拟次数。

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

s0 = 10. #刚买入的股价，初始状态
mu = 0.15
sigma = 0.2 
nSimu = 10000 # 模拟多少次
nStep = 244 # 每一次模拟走多少步, 1年只有244个交易日
dt = 1/nStep # 由于公式中的单位是年，而一年中有244个交易日，所以dt是理想化的每次股价变化间隔时间。

模拟10000次，得到一年后的股价分布。

price = np.empty((nSimu, nStep))
e = norm()
for i in range (nSimu):
    s = s0
    price[i, 0] = s
    for j in range(nStep):
        et = e.rvs() # 
        s = s + mu*s*dt + sigma*s*et*np.sqrt(dt) #股价的计算公式
        price[i, j] = s

画出一年后股价的分布。

plt.hist(price[:,-1], bins=30, density=True)
plt.ylabel('pdf')
plt.xlabel('Price')

5.1.2 采样

什么是采样？

采样（sampling）指的是从一个概率分布中抽取样本的过程。

采出来的样本对应的是概率分布的x轴
采出来很多样本之后，画histogram的图，应该近似概率分布

采样方法

一般来说，我们在python或matlab等语言中可以很轻松的借助内置函数来从一些常见分布中采样：

from scipy.stats import norm,gamma
a = norm.rvs() # 从高斯分布中采样
print(a)
b = gamma.rvs(1) # 从伽马分布中采样
print(b)

但关键问题是：假如不允许使用这些函数，我们又该如何采样呢？

接下来我们以高斯采样的例子介绍几种常见的采样方法。

问题：只给定均匀分布 $[0,1]$ 的样本，如何利用该均匀分布从一个高斯分布中采样?

方法1：逆变化采样（Inverse Transform Sampling）

算法：

从均匀分布 $[0,1]$ 采样得到 $u$
从高斯分布的累积分布函数的反函数 $x = F^{-1}(u)$ , $x$ 就是需要获得的样本

几乎所有的已知分布都可以用这种方法采样（matlab和python内部就是这么做的）。

但是，如果一个复杂分布，其累积分布函数的反函数未知，如何采样?

方法2：接受-拒绝采样（Accept-Reject Sampling）

采样前的准备：

先找一个提议分布（proposed distribution） $G$ ，概率密度函数 $g(x)$ 已知，可以是任意分布（例如高斯分布）
再找一个常数 $C$ ，使得对于任意的 $x$ ，都有 $C \cdot g(x) \geq p(x)$

算法：

从提议分布 $G$ 中进行采样，获取采样样本 $x$
从 $[0, 1]$ 的均匀分布中进行采样一个数值 $u$
如果 $u \leq \frac{p(x)}{C \cdot g(x)}$ ,那么我们就接受这个样本 $x$ ，保存下来；如果不满足，就拒绝并丢弃这个采样值，重来

以下是代码示例：

首先，设置目标采样分布的概率密度函数：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm,uniform
# 目标采样分布的概率密度函数(已知但是很复杂)
def p(x):
    return ((0.3 * np.exp(-(x-0.3)**2))+0.7*np.exp(-(x-2.)**2/0.3))/1.2113

选定 $C$ 值和提议分布函数 $g(x)$ ，并画出分布。

# 设定C值
C = 2.5

# g(x)为均值=1.4, 标准差=1.2的高斯分布

# 画出两个分布
x = np.arange(-4., 6., 0.01)
plt.plot(x, p(x), color='r', label='p(x)')
plt.plot(x, C*norm.pdf(x, loc=1.4, scale=1.2), color='b', label='C*g(x)')
plt.legend()
plt.xlabel('x')
plt.ylabel('pdf')

现在我们来采样:

sample = []
    
nSample = 10000
for i in range(nSample):
    x = np.random.normal(loc=1.4, scale=1.2) 
    u = np.random.rand()
    if u <= p(x)/(C*norm.pdf(x, loc=1.4, scale=1.2)):
        sample.append(x)

完成之后画出分布：

x = np.arange(-3., 5., 0.01)
plt.plot(x, p(x), color='r', label='p(x)')
plt.plot(x, C*norm.pdf(x, loc=1.4, scale=1.2), color='b', label='C*g(x)')
plt.hist(sample, color='k', bins=150, density=True)
plt.legend()
plt.ylabel('p(x)')

以上采样过程还可以写成矩阵形式

nSample = 10000
# norm和uniform函数都是产生一个正态或均匀分布的对象，可以搭配rvs函数等使用
uniform_rv = uniform(loc=0, scale=1) 
norm_rv = norm(loc=1.4, scale=1.2)
yy = norm_rv.rvs(size=nSample)
uu = uniform_rv.rvs(size=nSample)
mask = (uu <= p(yy)/(C*norm_rv.pdf(yy)))
sample = yy[mask]

接受拒绝采样的直观理解：

在每一个 $x$ 的点上, 都有 $p(y) \leq C \cdot g(y)$ ，即 $0 < \frac{f(y)}{C \cdot g(y)} \leq 1$
从图上看，就是蓝线和红线的比值。

但是，接受拒绝采样需要找到一个准确的常数 $C$ ，如果这个常数 $C$ 无法直接得到怎么办？

方法3：Metropolis-Hasting采样

算法：

从任意一个 $x_0$ 开始
根据当前的样本 $x_t$ 构建一个提议分布 $g(x | x_t)$ （例如高斯分布），然后采样得到 $x^*$
然后计算接受率 $\alpha = \min \left[ 1, \frac{p(x^*) \cdot g(x_t \mid x^*)}{p(x_t) \cdot g(x^* \mid x_t)} \right]$
从 $[0,1]$ 采样一个值 $u$ ,
- 如果 $u< \alpha$ , 接受该样本 $x_{t+1} = x^*$
- 否则，继续维持原来的样本 $x_{t+1} = x_t$
重复2-4步直到达到设定的样本总数

在MCMC的证明中会详细介绍该算法并解释它为什么成立，在这里，我们先来看看代码如何实现。

首先，设置目标分布的概率密度函数：

from scipy.stats import norm
import numpy as np 

# 目标采样分布的概率密度函数
def p(x):
    return ((0.3 * np.exp(-(x-0.3)**2))+0.7*np.exp(-(x-2.)**2/0.3))/1.2113
    
nSample = 100000
sample2 = []

y = 1 # 设置初始值

#进行MH采样
for i in range(nSample):
    
    y2 = norm(loc=y).rvs() # 从提议分布采样一个
    
    a = np.min([1, p(y2)/p(y)]) # 因为提议分布高斯是对称的, 所以可以上下约掉

    u = np.random.rand() # 从[0,1]均匀分布获取一个值
    
    if u <= a:
        sample2.append(y2) # 接受并跳转到新的sample
        y = y2
    else:
        sample2.append(y) # 保留原来那个sample不变
        
# 绘制采样结果
x = np.arange(-3., 5., 0.01)
plt.plot(x, p(x), color='r', label='p(x)')
a,_,_=plt.hist(sample2, color='k', bins=150, density=True)
plt.legend()
plt.xlabel('X')
plt.ylabel('p(x)')

5.1.3 马尔可夫链

什么是马尔可夫链？

未来只依赖于现在，而不依赖于过去。

马尔可夫链是状态空间中从一个状态到另一个状态的转换的随机过程，用来描述一系列状态之间的转移概率。在转移过程中，下一个状态 $S_{t+1}$ 只和当前状态 $S_t$ 有关，和再之前的状态 $S_{t-1}$ 无关。

P(x_{t+1}|x_t,x_{t-1},x_{t-2},...,x_t) = P(x_{t+1}|x_t)

马尔可夫链：状态转移、稳态分布和细致平稳条件

在马尔可夫链中，有两个核心概念：状态转移与稳态分布。

我们下面用一个情绪变化的例子来理解这两个概念。

假定人的情绪有三种状态：高兴、一般和悲伤。在每个时间点，人的情绪可以在这三种状态之间转换，且此刻的情绪只和上一个时刻的情绪状态有关，与再之前的情绪无关。

我们可以用状态转移矩阵 $Q$ 来描述这种转换。

在马尔可夫链中，随着状态转移过程，状态分布会发生什么样的变化呢？

假设在 $t_0$ 时刻的初始情绪概率分布为 $p_0$ ，状态转移矩阵为 $Q$ 。我们可以借助代码一步一步计算，来跟踪情绪的演化。

首先，我们定义一个情绪状态的初始分布 $p_0=[0.5,0.3,0.2]$

import numpy as np
# 定义初始状态p0
e = np.array([0.5, 0.3, 0.2])

然后，我们定义一个 $3\times 3$ 的转移矩阵 $Q$ ，第 $i$ 行代表从第 $i$ 种情绪状态转移到三种情绪状态的概率

import numpy as np
# 定义状态转移矩阵Q
Q = np.array([0.9, 0.075, 0.025, 0.15, 0.8, 0.05, 0.25, 0.25, 0.5])
Q = Q.reshape(3,3) # 得到3*3矩阵

接着，从初始状态开始，以转移矩阵进行40次状态转移：

# 状态转移
for i in range(40): 
    e = e.dot(Q)
    print(f'interation i={i}, {e[0]:.3f} {e[1]:.3f} {e[2]:.3f}\n')

我们可以得到如下结果：

interation i=0, 0.545 0.328 0.128
interation i=1, 0.572 0.335 0.094
interation i=2, 0.588 0.334 0.078
interation i=3, 0.599 0.331 0.070
interation i=4, 0.606 0.327 0.067
...
interation i=36, 0.625 0.313 0.063
interation i=37, 0.625 0.313 0.063
interation i=38, 0.625 0.313 0.063
interation i=39, 0.625 0.313 0.063

你发现了吗？经过多次状态转移，情绪状态的概率分布最后会保持不变： $p=[0.625,0.313,0.063]$ 。

那如果我们换一个初始条件呢？

# 定义初始状态
e = np.array([0.1, 0.1, 0.8])
# 状态转移
for i in range(40): 
    e = e.dot(Q)
    print(f'interation i={i}, {e[0]:.3f} {e[1]:.3f} {e[2]:.3f}\n')

输出：

interation i=0, 0.305 0.288 0.408
interation i=1, 0.420 0.355 0.226
interation i=2, 0.487 0.372 0.141
interation i=3, 0.530 0.369 0.101
interation i=4, 0.557 0.360 0.082
...
interation i=36, 0.625 0.313 0.063
interation i=37, 0.625 0.313 0.063
interation i=38, 0.625 0.313 0.063
interation i=39, 0.625 0.313 0.063

改变了初始条件后，经过状态转移过程，不仅分布依然稳定下来，而且得到了一模一样的分布！

由此，我们引出马尔可夫链的两个性质：

对于初始的状态分布，如果固定了状态转移矩阵（离散状态）或状态转移函数（连续状态），那么多次迭代之后会收敛到一个最终的稳定的分布。
无论初始状态分布如何，都会收敛到那个稳态分布。

也就是说，如果我们想要得到一个目标分布，只需要找到对应的状态转移矩阵，然后随机初始状态，通过不断迭代就能够达到收敛。

因此，我们最后介绍马尔可夫链的细致平稳条件：

满足以下条件的状态分布 $p$ 就是状态转移矩阵（或函数） $Q$ 的细致平稳分布：

p(x_1)*Q(x_2|x_1)=p(x_2)*Q(x_1|x_2)

马尔可夫链在MCMC中的作用

在MCMC中，我们旨在构造一个马尔可夫链，其稳态分布就是我们想要采样的目标后验分布。

具体来说，沿着马尔可夫链采样：

每一步只需要当前状态，不需要追踪过去整个状态历史。
无论目标分布是否已知，只要找到正确的转移矩阵（或函数），理论上就能收敛到目标分布。
细致平稳条件确保了我们设计的转移矩阵最终能够收敛到目标分布上。
可以在复杂、高维状态空间中缓慢但稳定地探索，避免维度灾难。

下面我们还是用前面提到的情绪变化的例子，示例如何用代码进行马尔可夫蒙特卡洛采样。

定义初始情绪状态：

# 定义一个初始状态
s = 0 # 3种情绪状态, s=0,1,2

设置烧录（burn-in）和采样的样本量（烧录的目的是保证取样时分布已经稳定，所以一般会取比较大的数）：

burnin = 1000000 # 被丢弃的burn-in samples
nSample = 1000000
samples = []

进行状态转移：

for i in range(burnin + nSample):
    u = np.random.rand()
    match s:
        case 0:
            if u <= 0.9:
                s = 0
            elif (0.9<u) & (u <= 0.9+0.075):
                s = 1
            elif 0.9+0.075 < u:
                s = 2
        case 1:
            if u <=0.15:
                s = 0
            elif (0.15<u) & (u <= 0.15+0.8):
                s = 1
            elif 0.15+0.8 < u:
                s = 2
        case 2:
            if u <=0.25:
                s = 0
            elif (0.25<u) & (u <= 0.25+0.5):
                s = 1
            elif 0.25+0.5 < u:
                s = 2
    samples.append(s)

得到采样分布：

bins = [ -0.5, 0.5, 1.5, 2.5 ]
a, _, _ = plt.hist(samples[burnin:], bins=bins)
plt.xlabel('State')
plt.ylabel('Count')
print(a/nSample)

采样得到的分布接近马尔可夫链的稳态分布[0.625, 0.313, 0.062]：

[0.617709 0.339236 0.043055]

5.1.4 MCMC：MH采样算法证明

到这里，我们已经把马尔可夫链蒙特卡洛采样拆成三部分分别解释了一遍。

马尔可夫链蒙特卡洛采样（Markov Chain Monte Carlo Sampling，MCMC）关键词含义：

马尔可夫链（Markov Chain）: 采样是沿着一条马尔可夫链进行状态转移
蒙特卡洛采样（Monte Carlo Sampling）: 用采样的方法来逼近某个分布

现在，我们来把前面讲的三部分结合起来。

还记得前面的Metropolis-Hasting采样吗？我们现在可以来解释为什么这样的算法是可行的了。

再回顾一下MH算法：

从任意一个 $x_0$ 开始
根据当前的样本 $x_t$ 构建一个提议分布 $g(x | x_t)$ (e.g. $g(x | x_t)$ 可以是一个高斯分布)，然后采样得到 $x^*$
然后计算接受率 $\alpha = \min \left[ 1, \frac{p(x^*) \cdot g(x_t \mid x^*)}{p(x_t) \cdot g(x^* \mid x_t)} \right]$
从 $[0,1]$ 采样一个值 $u$ ,
- 如果 $u< \alpha$ , 接受该样本 $x_{t+1} = x^*$
- 否则，继续维持原来的样本 $x_{t+1} = x_t$
重复2-4步直到达到预订的样本数量

这里有两个疑问：

提议分布 $g(x | x_t)$ 是由我们自己随意定的，凭什么随意定一个 $g$ 都可以能从 $p(x)$ 进行采样？
为什么接受率设置成这个形式之后，就可以对目标分布 $p(x)$ 进行采样？

从马尔可夫链的性质，我们知道，要从复杂的 $p(x)$ 中采样，就要找到与 $p(x)$ 对应的状态转移函数。怎么找呢？随便假定一个转移函数 $Q$ 显然是不行的，它不满足细致平稳条件，即：

p(x_i)*Q(x_j | x_i) \not= p(x_j)*Q(x_i | x_j)

所以我们要对以上的公式做一定的改造，让它满足细致平稳条件。

首先，引入 $a(x_j | x_i)$ , 使得上面的式子可以取等号, 即：

p(x_i)*Q(x_j | x_i)*a(x_j | x_i) = p(x_j)*Q(x_i | x_j)*a(x_i | x_j)

问题是什么样子的 $a(x_j | x_i)$ 可以满足该要求呢？其实很简单，只要满足下面条件：

a(x_j | x_i) = min(1, \frac{p(x_j)*Q(x_i | x_j)}{p(x_i)*Q(x_j | x_i)})

把 $a(x_j | x_i)$ 带入公式，

\begin{align*} p(x_i)*Q(x_j | x_i)*a(x_j | x_i) &= p(x_i)*Q(x_j | x_i)*min(1, \frac{p(x_j)*Q(x_i | x_j)}{p(x_i)*Q(x_j | x_i)}) \\ &= min(p(x_i)*Q(x_j | x_i), p(x_j)*Q(x_i | x_j)) \\ &= p(x_j)*Q(x_i | x_j) * min(\frac{p(x_i)*Q(x_j | x_i)}{p(x_j)*Q(x_i | x_j)}, 1) \\ &= p(x_j)*Q(x_i | x_j)*a(x_i | x_j) \end{align*}

如果我们把新的转移函数设置为 $M(x_j | x_i)=Q(x_j | x_i)*a(x_j | x_i)$ ，自然就满足细致平稳条件：

p(x_i)*M(x_j | x_i)= p(x_j)*M(x_i | x_j)

注意，这里的 $Q(x_j | x_i)$ 可以是任意的转移函数（例如高斯函数），以上式子恒成立。

回答前面两个疑问：

任意提议分布 $g(x | x_t)$ 就任意是转移函数 $Q(x_j | x_i)$ ，它的确无法直接满足细致平稳条件。
但是，任意提议分布 $g(x | x_t)$ 和接受率 $\alpha$ 结合起来形成了新转移函数 $M(x_j | x_i)$ ：先从任意提议分布中采样，然后接受率决定是否采纳新样本——新的转移函数满足细致平稳条件。

我们知道，认知心理学中概率模型的本质就是求解后验概率分布 $p(\theta | data)$ ，

p(\theta | data) = \frac{p(data | \theta) \cdot p(\theta)}{p(data)}

换句话来说，我们的目标分布是 $p(\theta|data)$ 。为了和前面的公式推导符号一致，我们把目标分布写为 $p(x|data)$ ，并带入公式，得到：

\begin{align*}a(x_j | x_i) &= min(1, \frac{p(x_j | data) \cdot Q(x_i | x_j)}{p(x_i | data) \cdot Q(x_j | x_i)}) \\ &= min(1, \frac{p(data | x_j) \cdot p(x_j)/p(data) \cdot Q(x_i | x_j)}{p(data | x_i) \cdot p(x_i)/p(data) \cdot Q(x_j | x_i)})\end{align*}

假设先验分布为均匀分布，则 $p(x_i) = p(x_j)$ ；假设提议分布是对称的（例如高斯分布，实际上前面说过的任意提议分布都可以），则 $Q(x_j | x_i) = Q(x_i | x_j)$ 。上式进一步简化为：

a(x_j | x_i) = min(1, \frac{p(data | x_j)}{p(data | x_i)})

可以看到，等式右边直接就是似然函数的比值。

如果想要进一步理解数学推导，可以参考以下资料：

博客资料: https://www.cnblogs.com/pinard/p/6625739.html
B站视频: https://www.bilibili.com/video/av32430563/

5.1.5 MCMC在心理学中的应用

MCMC在参数估计中的应用

在第四章的心理物理曲线部分，我们利用最大似然估计的方法拟合了心理物理曲线。这种方法的问题在于，如果中间的隐变量无法通过积分或求和消除，就很难求解。

我们通过手写MH算法来估计心理物理曲线。

实验任务是一个随机点运动任务。实验设置的一致性水平为[0.032, 0.064, 0.128, 0.256, 0.512]。已知该被试的行为反应数据，包含 $5个一致性水平 \times 1000个试次 = 5000 个试次$ 。我们的目标是用Weilbull函数拟合一位被试的心理物理曲线，并估计出82%阈值。

Weibull函数为

y=1-(1-\gamma)exp^{-(\frac{k*x}{\alpha})^\beta}

其中，

k=-log(\frac{1-g}{1-\gamma})^\frac{1}{\beta}

$\alpha$ 为该正确率下的阈值， $\beta$ 为函数的斜率， $\gamma$ 为基线概率水平的概率值， $g$ 为阈值对应的正确率。

该被试的心理物理曲线参数组合 $(\alpha,\beta,\gamma,g)$ 的真值为（0.25，3，0.5，0.82），其产生的行为数据如下。具体的数据模拟代码见第四章的心理物理曲线部分

首先，与最大似然估计类似，我们先写出关于 $\alpha$ 的负对数似然函数：

# 定义weibull函数
def weibullfun(x, alpha, beta, gamma, g):
    k = (-np.log((1-g)/(1-gamma)))**(1/beta)
    y = 1 - (1 - gamma)*np.exp(-(k*x/alpha) ** beta) # hack here to avoid the complex number of 
    return y

# 定义负对数似然函数
eps = np.finfo(float).eps #大于0的极小值，防止后续计算对数时出现0的情况
def negloglikeli(alpha, coh, data):
    gamma = 0.5
    g = 0.82
    beta = 3.

    prob = np.empty(coh.size)
    for i in range(coh.size):
        p = weibullfun(coh[i], alpha, beta, gamma, g)
        prob[i] = p*data[i]+(1-p)*(1-data[i])
    
    prob = prob*0.999+eps # 避免prob里面有0的情况, log(0)=-infinity
    return -np.log(prob).sum() 
    #prob取值在（0，1]，经过对数函数后会变为负数，取负的对数似然值能保证返回值是一个正数

使用MH算法对后验分布 $p(\alpha|data)$ 进行采样，选择对称的高斯分布进行提议分布，便可以使用公式 $a(x_j|x_i)=min(1,\frac{p(data|x_j)}{p(data|x_i)})$ 计算接受率。考虑到在代码实现时，直接计算似然的比值有些麻烦，可以两边取对数，改为计算两次采样参数的对数似然的差值：

log[a(x_j|x_i)]=min(0,log[p(data|x_j)]-log[p(data|x_i)])

MCMC采样从初始值0（可任意选取）开始，共采样2000次，剔除前1000次未平稳收敛的采样，保留后1000个有效样本。

from scipy.stats import norm

nSample = 1000
nBurnin = 1000
nTotalSample = nSample + nBurnin
stepsize = 1

samples = []
oldsample = 0

由于提议分布是在实数范围内进行采样，但需要估计的阈限 $\alpha$ 参数范围在(0,1)，我们需要进行一个尺度转换，通过logistic函数可以将实数范围投射到(0,1)之中。

trans_oldsample = 1/(1+np.exp(-oldsample))

每次采样时，会在均值为上一采样点的高斯分布中采样，进行尺度转换后通过负对数似然函数计算接受概率，并根据随机数决定是否接受本次采样。

for i in range(nTotalSample):
    newsample = norm.rvs(loc=oldsample,scale = stepsize)
    # 我们的参数范围在(0,1),但是我们转化到纯实数范围内采集
    trans_newsample = 1/(1+np.exp(-newsample))

    # 计算负对数似然之差
    delta = negloglikeli(trans_oldsample, cohTrial, data)-negloglikeli(trans_newsample, cohTrial, data)
    # 接受概率
    a = np.min([0,delta]) 
    #从[0,1]均匀分布随机获取一个值
    u = np.random.rand()
    
    if log(u) <= a: #接受并跳转到新的采样点
        samples.append(trans_newsample)
        trans_oldsample = trans_newsample
        oldsample = newsample
    else: #保留原来的采样点不变
        samples.append(trans_oldsample)

完成之后可以从样本计算得到所估计参数阈限 $\alpha$ 的均值，并画出有效采样的直方图。由于模拟数据与采样具有随机性，运行的结果可能有所不同。

print(np.mean(samples[1000:]))
plt.hist(samples[1000:], bins=20,density=True,range=[0.2,0.3])

0.24384956185678375

可以看到使用MH算法估计的参数 $\alpha$ 的均值，接近参数设置的真值0.25。

在采样数量较小时，分布可能很不稳定，可以通过提高采样数量来获得更高精度，但是运行速度会更加慢。

MCMC工具包选择指南

MCMC的算法本身和具体的模型无关，只是一种采样的方法，所以现在有许多工具包可以实现MCMC的算法。

目前的主流MCMC算法包括:

Metropolis-Hasting sampler
Gibbs sampler
Hamiltonian sampler
No-U-turn sampler

目前的主流Bayesian编程工具包包括：

Stan (基于C++, 有pystan, rstan接口，推荐)
PYMC5 (基于python，不方便写for循环)
Numpyro (基于python，不方便写for循环)
Jags, winbugs (已被淘汰，不推荐)

5.1.6 MCMC方法的优势与劣势

MCMC方法的优势

估计一个参数的分布而不是点估计，可以知道估计的误差
一旦生成模型很复杂，likelihood函数我们是无法直接通过解析解写出来的，那么 MCMC就成了唯一的方法

MCMC方法的劣势

因为采样必须要按照马尔可夫链依次序进行，没法进行平行加速，所以进程很慢
虽然最后得到了很多后验分布的样本，依然不太好做统计检验。

因此，在实际操作，尤其是对于计算效率要求比较高的环境中（比如深度学习模型训练和大规模贝叶斯推断），更多地会采用下节要介绍的变分推断等方法。

最后更新于2个月前