马尔科夫链蒙特卡洛采样

最后更新于6个月前

图片+应用+贝叶斯脑理论

写在前面的话：

事情要从这样一篇“不知所云”的文章讲起。（图片图片图片图片图片

虽然研究对象是熟悉的决策风格，也挂上了神经生物学的名头，并且成功地推送到了某心理学学生面前。但是显然，这些来自于“computational”的专业名词，诸如“hierarchy (hyper-group parameters)”“uninformed priors (i.e, uniform distributions)”“the Markov Chain Monte Carlo (MCMC) technique”“initial burn-in sequence”等等，对于心理学学生来说实在是如同天书。即使努力询问或是自学概率统计和机器学习，也难以打通学科间的壁垒，融会贯通到心理学上去。

困难虽有，不必着急，本章内容也许能帮助你迈出理解第一步。

在中我们已经对和有了详尽介绍，接下来要介绍的马尔科夫链蒙特卡洛采样（Markov Chain Mento Carlo Sampling, MCMC）和与后者同属于贝叶斯学派，但又有所不同，可以根据不同的标准进行划分。

不同的分类标准（树状图树状图树状图

按照参数是精确估计还是近似估计可以区分为精确推断（exact inference）和近似推断（approximate inference）两类。其中最大似然估计和最大后验估计属于前者，需要精确计算参数的概率分布，这类方法可以用于较为简单的模型；马尔科夫链蒙特卡洛采样和变分推断则属于近似推断，在参数不可计算或是计算成本过高时，使用近似方法推断参数的后验概率或统计特性。
按照参数估计是针对似然方程还是后验分布也可以区分。最大似然估计优化似然方程参数得到点估计，即单一估计值；最大后验估计优化后验参数也得到点估计，马尔科夫链蒙特卡洛采样和变分贝叶斯方法则获得完整后验分布（有可能以不同的形式）。

在认知神经科学中，最常用的是MLE, MAP和MCMC

相较MLE/MAP, MCMC的优势

MLE基于似然函数，很多时候似然函数里面包含多个隐变量和多重积分，无法得到解析形式，无法优化求解。MCMC不需要解析求解任何积分。
MLE只是关于参数的点估计。MCMC得到整个后验概率分布，可以得到参数估计的uncertainty
MLE很容易会收敛到局部最小值（local minima），MCMC一般不会。

相较MLE/MAP, MCMC的劣势

收敛速度慢
得到后验概率之后，不方便做统计以及其他分析。(如何根据sampling近似的后验概率做各种分析是一个热门话题！)

对于非machine learning专业的人来说，理解MCMC确实是非常大的挑战，但相应的，回报也很丰厚，能够理解MCMC可以说是理解probabilistic Bayesian modeling的一个里程碑。

理解MCMC的三个层次

完全理解所有数学推导和应用代码实现(excellent!)
部分理解背后数学推导但是可以实现采样算法(good!)
部分或者不理解背后数学推导但是可以借助例如stan软件实现采样算法(still good!)

学习认知神经科学的人基本需要做到2和3就可以了。

在学习过程中，不少心理学学生由衷地表示：“啊，数学好复杂，说人话”“好吧，好不容易看懂数学，和我(心理学)有什么关系”“似乎有点关系哦，那代码怎么整”。这些疑问在看到学科之间相融的巧妙时会自然消失，知识上的不足也会在学习中补齐，话不多说，我们开始正式学习马尔科夫链蒙特卡洛采样。

马尔科夫链

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

举个例子:

eg.1 假定人的情绪有三种状态：高兴、一般和悲伤。在某个时间点，人的情绪可以在这三种状态之间转换。我们可以用状态转移矩阵 $P$ 来描述这种转换。假设在 $t_0$ 时刻的初始情绪分布为 $p_0=[0.5,0.3,0.2]$ ，我们可以通过一步一步的计算来跟踪情绪的演化。（图片图片图片

首先定义转移矩阵 $P$ :

import numpy as np
# define transition matrix
P = np.array([0.9, 0.075, 0.025, 0.15, 0.8, 0.05, 0.25, 0.25, 0.5])
P = P.reshape(3,3) # 得到3*3矩阵 
print(P)

自定义初始分布并开始运行：

# initial states
e = np.array([0.5, 0.3, 0.2])

for i in range(40): 
    e = e.dot(P)
    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

可以看到情绪状态的分布最终稳定。

如果我们改变初始条件呢？

# initial states
e = np.array([0.1, 0.1, 0.8])

for i in range(40): #
    e = e.dot(P)
    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

改变了初始条件后，经过状态转移过程得到了相同的分布！

由此我们可以得出离散分布下的马尔科夫链性质：

对于一个初始的离散分布 $e_0$ ，如果固定了一个转移矩阵 $P$ ，那么多次迭代之后我们会发现收敛到一个最终的稳定的分布
无论初始状态分布如何，都会收敛到那个稳态分布

该性质可以推及连续分布，对比来看：
- 离散分布， $p(x)$ 是概率质量函数， $Q(x_2|x_1)$ 是一部分状态转移矩阵，代表从状态 $x_1$ 转移到状态 $x_2$ 的条件概率；
- 连续分布， $p(x)$ 是概率密度函数， $Q(x_2|x_1)$ 是一部分状态转移函数。

马尔科夫链的细致平稳条件

经推导（略），满足以下条件的分布 $p$ 就是转移矩阵 $Q$ 的细致平稳分布。

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

马尔科夫链中的采样

一个稳态分布对应一个转移矩阵，如果不利用原分布的概率质量函数，可以从任何一个状态开始，然后依照和稳态分布对应的转移矩阵或者转移函数进行状态转移,就可以达到从稳态分布采样的目的（迭代够多会得到稳定分布，无论初始状态）。

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

举个例子：

eg.2 就像eg.1 一样我们首先定义初始状态：

# 定义一个初始状态
t = 0 # 注意这里t是状态不是分布, t=0,1,2

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

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

进行状态转移。

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

可以直观地看出采样的分布。

a, _, _ = plt.hist(samples[burnin:], bins=3)
plt.ylabel('Count')
print(a/nSample)

应用：

如何从一个离散分布中采样呢？以 $[0.5, 0.3, 0.2]$ 为例，我们可以直接设置 $[0,0.5], [0.5,0.8], [0.8,1.0]$ 的区间，再利用random函数随机取样，代码不再写出。

那不知道概率质量函数的时候如何从中采样呢？

利用Markov的性质，如果我们找到其对应的转移矩阵，只要随意指定一个初始状态，然后沿着状态转移矩阵来采样转换状态，最终会达到稳态，也可以达到采样的目的。

蒙特卡洛方法

蒙特卡洛（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$ 会更加接近真实值

虽然无法直接计算圆的面积，但是可以在大量取点后，用 $\frac{m}{n}$ 来估计 $\pi$ 的值。

而我们经常要在一个分布上做积分，比如求期望，方差，做其他积分，但是一旦分布的概率函数比较复杂，不能做积分，就很难办。这时候就可以尝试绕过概率函数和分布函数，通过大量采样的方法直接得到统计量。

举个例子：

如果离散函数求期望或者连续函数求积分时，分布函数或者概率密度函数过于复杂，比如分布函数为 $P(x) = \frac{1}{1.2113} \left( 0.3 \cdot e^{-\frac{(x - 0.5)^2}{2 \cdot 0.1^2}} + 0.7 \cdot e^{-\frac{(x - 2.5)^2}{2 \cdot 0.5^2}} \right)$ ，概率密度函数为 $f(x) = e^{-0.5 \cdot x^2} + \sin(x) + x^3$ ，难以直接求解积分。

但是我们可以从 $p(x)$ 中取样 $x_1, x_2, x_3, ...x_i$ ，再利用大数定理（在本书中不做详细解释）求期望。

在离散分布中：

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

在连续分布中：

\mathbb{E}[X] = \int X_i\cdot p(x_i)\,d(x)\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}

原本需要求出股票处于某个价格的概率后求期望，运用蒙特卡洛方法则只需要模拟无数次后求均值。

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

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')

采样

在中其实已经提及到了什么是采样（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]$ 采样得到 $x$
从高斯分布的积累分布函数CDF的反函数 $y = f^{-1} CDF(x)$ , $y$ 就是需要获得的样本

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

可是如果一个复杂的分布，不知道其积累分布函数（CDF）的反函数，如何采样?

举个例子：给定一个均匀分布 $[0,1]$ ，如何从以下复杂分布中采样?

P(x) = \frac{1}{1.2113} \left( 0.3 \cdot e^{-\frac{(x - 0.5)^2}{2 \cdot 0.1^2}} + 0.7 \cdot e^{-\frac{(x - 2.5)^2}{2 \cdot 0.5^2}} \right)

方法2：接受-拒绝采样(Accept-Rejection sampling)

过程可以分为两步：

采样前的准备：

先找到一个建议分布(proposed distribution) $G$ 的概率密度函数是 $g(x)$ ，可以是任意分布(e.g. 高斯分布)
再找到一个常数 $C$ , 使得对于任意的 $x$ , 都有 $C \cdot g(x) \geq p(x)$

算法：

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

用代码来展示过程。

首先导入目标采样分布的概率密度函数。

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

可以画出该概率密度函数，如图。

x = np.arange(-4., 6., 0.01)
plt.plot(x, p(x), color='r', label='p(x)')
plt.legend()

选定 $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 = np.arange(-3., 5., 0.01)
plt.plot(x, p(x), color='r', label='p(x)')
plt.plot(x, C*norm_rv.pdf(x), color='b', label='C*g(x)')
plt.hist(sample, color='k', bins=150, density=True)
plt.legend()
plt.ylabel('p(x)')

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

在每一个 $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)$ (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步直到达到预订的样本数量

在中会用到这种算法并解释它为什么成立，我们先来看看代码如何实现。

首先表示出目标分布的概率密度函数。

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

进行MH采样。

nSample = 100000
sample2 = []

y = 1 # 设置初始值

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.ylabel('pdf(x)')

MCMC的证明

可以从字面意义上直接拆分MCMC，而这正是前几节我们做的事。

Markov Chain Monte Carlo Sampling 关键词的含义

Markov Chain: 我们采样是沿着一条Markov链在进行状态转移
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)$ 进行采样？

下面我们来进行证明。

前面知道, 根据Markov采样定理, 我们要 $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)})

把公式3带入公式2,

\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*} \tag{4}

可以发现，如果另一个新的转换函数 $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)$ 可以是任意的转换函数(e.g. 高斯函数)，上面的式子恒成立。

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

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

根据贝叶斯公式我们代进去：

\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/

小结（图片图片图片图片

MCMC sampling在心理学中的应用

作为参数估计的工具

（paper图片图片图片图片

回看开头的这篇文献，+术语解释

在认知计算中的应用

方法1：用MLE估计心理物理曲线的参数

详细代码见使用这种方法的缺点也很明显，中间的隐变量最好是能解析解积分消除，否则求解很麻烦。

方法2：手写MH算法来求解心理物理曲线的参数

像MLE方法中一样，我们画出韦伯曲线并收集数据后，写出负对数似然函数，假定threshold $\alpha$ 是需要估计的自由参数（这部分的代码省略，具体见）。

首先设置MCMC的采样数量，包括实际的采样数nSamples和丢弃的Burn-in阶段样本数nBurnin。stepsize是采样过程中步长的设定，它决定了每次采样新值与旧值的变化幅度。oldsample是MCMC采样的初始值，因为心理物理曲线中阈值在 $(0,1)$ 之间，因此要将初始值映射到 $(0,1)$ 的范围内。

rom numpy import log
from scipy.stats import norm

nSample = 1000
nBurnin = 1000
nTotalSample = nSample + nBurnin

stepsize = 1

samples = []

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

在每次迭代中，从正态分布中生成新的候选采样值newsample，并将其映射到 $(0,1)$ 。

for i in range(nTotalSample):
    newsample = norm.rvs(loc=oldsample, scale=stepsize)
    trans_newsample = 1 / (1 + np.exp(-newsample))

计算似然差异并决定是否接受新样本。

delta = negloglikeli(trans_oldsample, cohTrial, data) - negloglikeli(trans_newsample, cohTrial, data)

if delta > 0:  # 如果新样本比旧样本好，接受新样本
    samples.append(trans_newsample.copy())
    trans_oldsample = trans_newsample
    oldsample = newsample
else:
    tt = np.random.rand()
    if log(tt) < delta:  # 按一定概率接受差的样本
        samples.append(trans_newsample.copy())
        trans_oldsample = trans_newsample
        oldsample = newsample
    else:
        samples.append(trans_oldsample.copy())  # 保留旧样本

画出参数的分布。

plt.hist(samples[1000:])

（图片图片图片图片

可以求期望获得最终的 $\alpha$ 估计值。在bin较小时，分布可能很不稳定，可以通过提高采样数量来获得更高精度，但是运行速度会更加慢。

实现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 (已被淘汰，不推荐)

MCMC方法的优势与劣势

MCMC方法的优势

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

但MCMC方法的劣势也很明显，因为采样必须要按照马尔科夫链依次序进行，没法进行平行加速，所以进程很慢；此外，虽然最后得到了很多后验分布的样本，依然不太好做统计检验。

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

“贝叶斯脑”假说（The “Bayesian brain” Hypothesis）

首先定义转移矩阵 $P$ :

对于一个初始的离散分布 $e_0$ ，如果固定了一个转移矩阵 $P$ ，那么多次迭代之后我们会发现收敛到一个最终的稳定的分布

离散分布， $p(x)$ 是概率质量函数， $Q(x_2|x_1)$ 是一部分状态转移矩阵，代表从状态 $x_1$ 转移到状态 $x_2$ 的条件概率；

连续分布， $p(x)$ 是概率密度函数， $Q(x_2|x_1)$ 是一部分状态转移函数。

经推导（略），满足以下条件的分布 $p$ 就是转移矩阵 $Q$ 的细致平稳分布。

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

如何从一个离散分布中采样呢？以 $[0.5, 0.3, 0.2]$ 为例，我们可以直接设置 $[0,0.5], [0.5,0.8], [0.8,1.0]$ 的区间，再利用random函数随机取样，代码不再写出。

举个经典的例子来描述蒙特卡洛方法：估算圆周率 $\pi$

生成随机样本：在正方形内随机生成 $n$ 个点，记录落在圆内的点数 $m$

圆的面积与正方形面积的比值是 $\frac{\pi r^2}{(2r)^2} = \frac{\pi}{4}$

通过 $\frac{m}{n}$ 估计这个比值，得到 $\pi \approx 4 \times \frac{m}{n}$

增加样本数量：随着样本数量 $n$ ，估计值 $\pi$ 会更加接近真实值

虽然无法直接计算圆的面积，但是可以在大量取点后，用 $\frac{m}{n}$ 来估计 $\pi$ 的值。

但是我们可以从 $p(x)$ 中取样 $x_1, x_2, x_3, ...x_i$ ，再利用大数定理（在本书中不做详细解释）求期望。

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

\mathbb{E}[X] = \int X_i\cdot p(x_i)\,d(x)\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}

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

从均匀分布 $[0,1]$ 采样得到 $x$

从高斯分布的积累分布函数CDF的反函数 $y = f^{-1} CDF(x)$ , $y$ 就是需要获得的样本

举个例子：给定一个均匀分布 $[0,1]$ ，如何从以下复杂分布中采样?

P(x) = \frac{1}{1.2113} \left( 0.3 \cdot e^{-\frac{(x - 0.5)^2}{2 \cdot 0.1^2}} + 0.7 \cdot e^{-\frac{(x - 2.5)^2}{2 \cdot 0.5^2}} \right)

先找到一个建议分布(proposed distribution) $G$ 的概率密度函数是 $g(x)$ ，可以是任意分布(e.g. 高斯分布)

再找到一个常数 $C$ , 使得对于任意的 $x$ , 都有 $C \cdot g(x) \geq p(x)$

从建议分布 $G$ 中进行采样，获取采样样本 $y$

从 $[0, 1]$ 的均匀分布中进行采样，采取一个数值 $u$

如果 $u \leq \frac{p(y)}{C \cdot p(y)}$ ,那么我们就接受这个样本 $y$ ，保存下来; 如果不满足，就拒绝并丢弃这个采样值，重来。

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

在每一个 $x$ 的点上, 都有 $p(y) \leq C \cdot g(y)$ ，即 $0 < \frac{f(y)}{C \cdot g(y)} \leq 1$

可是接受拒绝采样需要找到一个准确的常数 $C$ , 如果这个常数 $C$ 不好直接得到怎么办?

从任意一个 $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$

从任意一个 $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$

这里的提议分布 $g(x | x_t)$ 是由我们自己随意定的，凭什么随意定一个 $g$ 都可以能从 $p(x)$ 进行采样？

为什么接受率定成这个形式之后就可以对 $p(x)$ 进行采样？

前面知道, 根据Markov采样定理, 我们要 $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)})

\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*} \tag{4}

可以发现，如果另一个新的转换函数 $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)$ 可以是任意的转换函数(e.g. 高斯函数)，上面的式子恒成立。

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

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

\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)})

像MLE方法中一样，我们画出韦伯曲线并收集数据后，写出负对数似然函数，假定threshold $\alpha$ 是需要估计的自由参数（这部分的代码省略，具体见）。

在每次迭代中，从正态分布中生成新的候选采样值newsample，并将其映射到 $(0,1)$ 。

可以求期望获得最终的 $\alpha$ 估计值。在bin较小时，分布可能很不稳定，可以通过提高采样数量来获得更高精度，但是运行速度会更加慢。