（终于把第十章读完了，这一章应该相对轻松。但这两天状态有待调整，所以没咋认真读）
目录
第十章介绍了确定性估计deterministic approximations，本章介绍蒙特卡洛方法。（关于蒙特卡洛方法和拉斯维加斯方法的区别，翻西瓜书~）
有些应用虽然关心非观测变量的后验分布，但大多数情况下，后验分布主要用于是求该分布下的期望
11.1 Basic Sampling Algorithms
这一节介绍一些从给定分布中生成随机变量的基本方法。并且假定算法生成 ( 0 , 1 ) (0, 1) (0,1)之间的均匀分布的伪随机数 z z z
P526 标准概率分布
使用某个函数
f
f
f对
z
z
z进行变换，即
y
=
f
(
z
)
y=f(z)
y=f(z)，
y
y
y上的概率分布为
选择一个函数
f
(
z
)
f(z)
f(z)使产生出的
y
y
y具有某种具体的分布形式
p
(
y
)
p(y)
p(y)，对上式积分得到
所以
y
=
h
−
1
(
z
)
y=h^{-1}(z)
y=h−1(z)
对于多个变量的推广，涉及到Jacobian行列式变换
例如Box-Muller方法生成二维高斯分布
先从均匀分布采样
z
1
,
z
2
∈
(
−
1
,
1
)
z_1, z_2\in (-1, 1)
z1,z2∈(−1,1)
然后丢弃
z
1
2
+
z
2
2
⩽
1
z_1^2 + z_2^2 \leqslant 1
z12+z22⩽1的点对，产生单位圆内的一个均匀分布，且
p
(
z
1
,
z
2
)
=
1
π
p(z_1, z_2)=\frac{1}{\pi}
p(z1,z2)=π1，对于每对
z
1
,
z
2
z_1, z_2
z1,z2，计算
其中
r
2
=
z
1
2
+
z
2
2
r^2=z_1^2 + z_2^2
r2=z12+z22. 这样得到的联合分布为
y
\bm y
y服从高斯分布，可以使用
μ
+
L
y
\bm \mu + \bm L \bm y
μ+Ly得到均值为
μ
\bm \mu
μ，协方差为
L
L
T
\bm L \bm L^T
LLT的高斯分布（如果要指定协方差为
Σ
\bm \Sigma
Σ，可以使用Cholesky分解，得到
Σ
=
L
L
T
\bm \Sigma = \bm L \bm L^T
Σ=LLT
P528 拒绝采样
我们假定
其中
p
~
(
z
)
\tilde p(z)
p~(z)好算，但是
Z
p
Z_p
Zp未知。为了执行拒绝采样，需要设计另一个提议分布proposal distribution. 引入一个常数
k
k
k，使得对所有的
z
z
z，都有
k
q
(
z
)
⩾
p
~
(
z
)
kq(z)\geqslant \tilde p(z)
kq(z)⩾p~(z)，函数
k
q
(
z
)
kq(z)
kq(z)叫做比较函数。如图所示
首先从
q
(
z
)
q(z)
q(z)中生成一个
z
0
z_0
z0，接着在区间
[
0
,
k
q
(
z
0
)
]
[0, kq(z_0)]
[0,kq(z0)]的均匀分布中生成
u
0
u_0
u0，这对于随机数在函数
k
q
(
z
)
kq(z)
kq(z)曲线下方是均匀分布。最后，如果
u
0
>
p
~
(
z
0
)
u_0> \tilde p(z_0)
u0>p~(z0)，那么样本被拒绝，否则
u
0
u_0
u0被保留。如上图阴影部分被拒绝。剩下的点在曲线
p
~
(
z
)
\tilde p(z)
p~(z)下方是均匀分布。因此
z
z
z服从
p
(
z
)
p(z)
p(z)
- 一个样本被接受的概率为
于是，常数 k k k应尽可能小，并处处不小于 p ~ ( z ) \tilde p(z) p~(z)
例如考虑从Gamma分布中采样
对均匀分布 y y y，使用 z = b tan y + c z=b\tan y+c z=btany+c进行变换，给出服从如下分布的随机数
最小拒绝率在如下条件得到， c = a − 1 , b 2 = 2 a − 1 c=a-1,b^2=2a-1 c=a−1,b2=2a−1，并将常数 k k k选的尽可能小，同时满足 k q ( z ) ⩾ p ~ ( z ) kq(z)\geqslant \tilde p(z) kq(z)⩾p~(z)，如下图所示
P530 可调节的拒绝采样Adaptive rejection sampling
许多情形中，构造合适的
q
(
z
)
q(z)
q(z)很困难。另一种确定函数形式的方法为基于
p
(
z
)
p(z)
p(z)直接构建界限函数envelope function. 对于
p
(
z
)
p(z)
p(z)是对数上凸函数，也即
ln
p
(
z
)
\ln p(z)
lnp(z)的导数非增时，构造方法如下图所示，算了几个格点
界限函数本身是由一个分段指数分布组成，形式为
如果一个样本被接受，那么就是采出来的一个样本，如果被拒绝，那么它会被并入格点集合中，算出一条新的切线，界限函数被优化，随着格点数量的增加，界限函数对所求概率分布的近似效果越来越好，拒绝的概率就会减小
- 该算法的一个变体，不用计算导数（Gilks,1992）
- 也可扩展到不是对数上凸函数的概率分布中，只需将每个拒绝采样的步骤中使用Metropolis-Hasting阶梯函数，从而产生adaptive rejection Metropolis sampling（Gilks et al.,1995）
- 高维空间下，接受率会迅速下降（一般是随维度指数级下降）。而且真实数据分布往往多峰，很难找到一个好的提议分布proposal distribution。书上用高维高斯分布举了个例子
- 对高维更复杂的算法来说，拒绝采样起着子过程的作用
P532 重要性采样
重要性采样则是在假定
p
(
z
)
p(\bm z)
p(z)采不到，但对于任意
z
\bm z
z都是可以算的前提下，算一个关于
f
(
z
)
f(\bm z)
f(z)的期望
与拒绝采样方式相同，引入容易采样的提议分布
q
(
z
)
q(\bm z)
q(z)
经常情况下，
p
(
z
)
=
p
~
(
z
)
/
Z
p
p(\bm z)=\tilde p(\bm z)/Z_p
p(z)=p~(z)/Zp中只有
p
~
(
z
)
\tilde p(\bm z)
p~(z)能计算，而
Z
p
Z_p
Zp未知。类似地，采用重要性采样
q
(
z
)
=
q
~
(
z
)
/
Z
q
q(\bm z)=\tilde q(\bm z)/Z_q
q(z)=q~(z)/Zq
其中，可以估计
从而
- p ( z ) p(\bm z) p(z)和 q ( z ) q(\bm z) q(z)越接近，该采样方法越成功
- 如果 p ( z ) f ( z ) p(\bm z)f(\bm z) p(z)f(z)变化剧烈，并集中在 z \bm z z的小区域，此时重要性权重 { r l } \{ r_l \} {rl}由几个具有较大值的权重控制，剩余权重较小。书上还对重要性采样的病态性有一些其他分析
- 如果是图模型，可以把观测值先设定到结点上，然后剩下的结点从均匀分布中采样，这样的到 q ~ ( z ) \tilde q(\bm z) q~(z)
- 似然加权采样likelihood weighted sampling，不写了，贴个公式
P534 采样-重要性-重采样 Sampling-importance-resampling
拒绝采样的常数
k
k
k并不好选，sampling-importance-resampling (SIR)方法则避免了选
k
k
k
该方法有两阶段，第一阶段，从
q
(
z
)
q(\bm z)
q(z)中采样
z
(
1
)
,
…
,
z
(
L
)
\bm z^{(1)}, \dots,\bm z^{(L)}
z(1),…,z(L)，第二阶段用式（11.23）构造
w
1
,
…
,
w
L
w_1, \dots, w_L
w1,…,wL. 最后，
L
L
L个关于
(
z
(
1
)
,
…
,
z
(
L
)
)
(\bm z^{(1)}, \dots,\bm z^{(L)})
(z(1),…,z(L))以概率
(
w
1
,
…
,
w
L
)
(w_1, \dots,w_L)
(w1,…,wL)进行采样。
- 当 L → ∞ L\to \infty L→∞时，分布变成了 p ( z ) p(\bm z) p(z)，书上给证了
- 当 q = p q=p q=p时， w n = 1 / L w_n=1/L wn=1/L
P536 采样和EM算法
采样方法可以近似EM算法中的E步，如果E步的解析解不好求的话。
在M步的优化目标是
可以通过采样近似
然后用通常方法该方法叫做Monte Carlo EM算法
一个特殊的形式是stochastic EM，在E步每次只采一个样本
P537 数据增广算法 data augmentation algorithm
假定从最大似然方式转换到贝叶斯方式，这里假设从
p
(
Z
∣
θ
,
X
)
p(\bm Z|\bm \theta, \bm X)
p(Z∣θ,X)和
p
(
θ
∣
Z
,
X
)
p(\bm \theta | \bm Z, \bm X)
p(θ∣Z,X)中采样相对简单
I(imputation step)P(posterior step) Algorithm
- I-step: 先从
p
(
Z
∣
X
)
p(\bm Z|\bm X)
p(Z∣X)采样
具体方法为从 p ( θ ∣ X ) p(\bm \theta|\bm X) p(θ∣X)采样 θ ( l ) \bm \theta^{(l)} θ(l)，然后从 p ( Z ∣ θ , X ) p(\bm Z|\bm\theta, \bm X) p(Z∣θ,X)采样 Z ( l ) \bm Z^{(l)} Z(l)，其中 l = 1 , … , L l=1,\dots,L l=1,…,L - P-step:
使用 { Z ( l ) } \{\bm Z^{(l)}\} {Z(l)}，估计后验分布 θ \bm \theta θ
这算法感觉是在交替估计两个边缘后验 Z \bm Z Z和 θ \bm \theta θ
11.2 Markov Chain Monte Carlo
上一节的rejection sampling和importance sampling，有很多限制，尤其是在高维
这一节的MCMC则是一个通用的框架，允许从更大类的分布中采样。并且能很好对应样本空间维度增长
MCMC中，提议分布为
q
(
z
∣
z
(
τ
)
)
q(\bm z|\bm z^{(\tau)})
q(z∣z(τ))，其中
z
(
τ
)
\bm z^{(\tau)}
z(τ)是当前状态的记录。这样能从马尔可夫链中采样
z
(
1
)
,
z
(
2
)
,
…
\bm z^{(1)},\bm z^{(2)},\dots
z(1),z(2),…，另一方面，假定
p
(
z
)
=
p
~
(
z
)
/
Z
p
p(\bm z)=\tilde p(\bm z)/Z_p
p(z)=p~(z)/Zp，其中
p
~
\tilde p
p~易算。算法的每一轮，从提议分布生成一个
z
∗
\bm z^*
z∗，然后以恰当的准则接受
P538 Metropolis算法
基础的Metropolis算法，假定了
q
(
z
A
∣
z
B
)
=
q
(
z
B
∣
z
A
)
q(\bm z_A|\bm z_B)=q(\bm z_B|\bm z_A)
q(zA∣zB)=q(zB∣zA). 其中接受率为
如果候选样本被接受，那么
z
(
τ
+
1
)
:
=
z
∗
\bm z^{(\tau+1)}:=\bm z^*
z(τ+1):=z∗；否则
z
(
τ
+
1
)
:
=
z
(
τ
)
\bm z^{(\tau + 1)}:=\bm z^{(\tau)}
z(τ+1):=z(τ)，这和拒绝采样不同
- 只要 q ( z A ∣ z B ) q(\bm z_A|\bm z_B) q(zA∣zB)是正的（充分不必要条件），当 τ → ∞ \tau \to \infty τ→∞时， z ( τ ) \bm z^{(\tau)} z(τ)的分布会倾向于 p ( z ) p(\bm z) p(z)，不过 z ( 1 ) , z ( 2 ) , … \bm z^{(1)}, \bm z^{(2)},\dots z(1),z(2),…不独立
考察一个随机游走的例子
如果
z
(
1
)
=
0
z^{(1)}=0
z(1)=0，可以求出
E
[
z
(
τ
)
]
=
0
\mathbb E[z^{(\tau)}]=0
E[z(τ)]=0，
E
[
(
z
(
τ
)
)
2
]
=
τ
/
2
\mathbb E[(z^{(\tau)})^2]=\tau / 2
E[(z(τ))2]=τ/2，所以走过的平均距离正比于
τ
\tau
τ的平方根，这是随机游走的典型性质，这说明随机游走的搜索低效。设计MCMC的一个中心目标是避免随机游走
P539 马尔可夫链
一阶马尔可夫链
转移概率为
T
m
(
z
(
m
)
,
z
(
m
+
1
)
)
=
p
(
z
(
m
+
1
)
∣
z
(
m
)
)
T_m(\bm z^{(m)}, \bm z^{(m+1)})=p(\bm z^{(m+1)}|\bm z^{(m)})
Tm(z(m),z(m+1))=p(z(m+1)∣z(m))
如果转移概率对所有
m
m
m都相同，则称之为同质的homogeneous
当满足
概率分布不变。一个给定的马尔可夫链可能有多个不变的概率分布
易证，当满足特定概率分布的细节平衡detailed balance的概率转移，会使概率分布具有不变性
满足细节平衡的马尔可夫链被称为可翻转的reversible
所以，我们要构造马尔可夫链，并从一个给定分布采样。实际中，经常构造一组基
B
1
,
…
,
B
K
B_1, \dots, B_K
B1,…,BK
并对混合系数满足
α
1
,
…
,
α
K
\alpha_1, \dots, \alpha_K
α1,…,αK满足
α
k
⩾
0
\alpha_k \geqslant 0
αk⩾0，且
∑
k
α
k
=
1
\sum_k \alpha_k=1
∑kαk=1. 另一种构造方法为
P541 Metropolis-Hastings算法
从
q
k
(
z
∣
z
(
τ
)
)
q_k(\bm z|\bm z^{(\tau)})
qk(z∣z(τ))采样
z
∗
\bm z^*
z∗，然后用如下概率接受
这里，
k
k
k标记出可能的转移集合中的成员。对于对称的提议分布，会退化为式11.33的标准Metropolis算法
可以证明
p
(
z
)
p(\bm z)
p(z)对于Metropolis-Hastings算法定义的马尔可夫链是一个不变的概率分布
- 提议分布经常选择以当前状态为中心的高斯，这里有个trade-off，如果方差太小，接受率高，但是随机游走得慢；如果方差过大，拒绝率很高，因为会到达
p
(
z
)
p(\bm z)
p(z)很低的位置。如图所示
11.3 Gibbs Sampling
Gibbs采样是一种简单通用的MCMC算法，并且是Metropolis-Hastings的特例
对于
p
(
z
)
=
p
(
z
1
,
…
,
z
M
)
p(\bm z)=p(z_1, \dots, z_M)
p(z)=p(z1,…,zM)，每次（随机或轮流）选一个进行更新
Gibbs采样中
- p ( z ) p(\bm z) p(z)是invariant的
- 该采样方法要满足各态历经性ergodic，即马尔可夫链最终收敛的分布与初始状态无关。最终收敛的状态叫做均衡分布equilibrium。各态经历性的一个充分条件是没有条件概率处处为0，也即 z \bm z z空间中任何一点都可以从其他一点有限步达到（习题11.12是一个例子）。如果不满足，则必须显式证明ergodic。
下面说明Gibbs采样是一种特殊的Metropolis-Hastings算法
如果看作
q
k
(
z
∗
∣
z
)
=
p
(
z
k
∗
∣
z
∖
k
)
q_k(\bm z^*|\bm z)=p(z_k^*|\bm z_{\setminus k})
qk(z∗∣z)=p(zk∗∣z∖k)，则
采样的例子如图所示
- 在该例子中，可以通过旋转坐标系，获得不相关的两个变量，但是实际应用中，很难这么搞
- Gibbs采样中一种减少随机游走行为的方式叫做over-relaxation过松弛，该方法应用于条件分布是高斯分布的情况（联合分布可以不是高斯）
在Gibbs采样中的每一步里， z i z_i zi的条件分布有均值 μ i \mu_i μi和方差 σ i 2 \sigma_i^2 σi2，而 z i z_i zi被代替为
其中 ν \nu ν是一个标准高斯的随机变量， α \alpha α则介于-1到1之间，当 α = 0 \alpha=0 α=0时，上式相当于标准吉布斯采样。当 α < 0 \alpha < 0 α<0时，该方法偏向均值的相反方向 - 在概率分布能表示成图的情况下，各个节点的条件概率只依赖于对应马尔可夫毯中的变量
- 如果Gibbs采样中每一步，不是采样，而是选最大似然的点估计，那么得到迭代条件峰值算法iterated conditional models(ICM) algorithm. 这可以看作是Gibbs采样的贪心近似
- 如果一次不是改变一个变量，而是更新一组，则得到了分块blocking Gibbs采样算法
11.4 Slice Sampling
Metropolis算法中对步长十分敏感。切片采样slice sampling提供了一种可调节步长。该方法需要能够计算未归一化分布
p
~
(
z
)
\tilde p(\bm z)
p~(z)
考虑一元变量的情况，切片采样需要使用
u
u
u对
z
z
z进行增广，然后从联合分布中采样。目标是从下列分布中采样
（这感觉像是从概率密度函数的二维集合区域中采样）
u
u
u和
z
z
z交替采样。主要是固定
u
u
u的时候怎么采
z
z
z，这里要设计一个包含
z
(
τ
)
z^{(\tau)}
z(τ)的动态长度的区间。一开始给长一些，从中采样，如果被拒绝了，就从拒绝点那里设为边界点，缩小，再采样。
具体细节翻书吧
- 切片采样可以用于多元分布，方法是按照吉布斯采样的方式重复对每个变量采样，这要求对元素 z i z_i zi，能够计算正比于 p ( z i ∣ z ∖ i ) p(z_i|\bm z_{\setminus i}) p(zi∣z∖i)的函数
11.5 The Hybrid Monte Carlo Algorithm
Metropolis算法有一个弊端，就是随机游走行为
P548 动力系统Dynamical systems
该方法源自哈密顿动力学Hamiltonian dynamics下物理系统的演化模拟. 位置为
z
\bm z
z，动量变量（其实就是速度吧）为
r
\bm r
r，两者形成相空间phase space
(
z
,
r
)
(\bm z, \bm r)
(z,r)
采样时，把连续的过程离散化，对于一个完整的时间间隔
τ
\tau
τ，用
ϵ
\epsilon
ϵ作为一个步长进行离散，这种近似误差会在
ϵ
→
0
\epsilon \to 0
ϵ→0时消失。具体更新为（蛙跳算法）
P552 混合Monte Carlo
混合Monte Carlo结合哈密顿动力系统和Metropolis算法
对动量变量
r
\bm r
r随机更新，和哈密顿动力系统的更新（蛙跳算法），两步交替完成。在每次应用蛙跳算法后，基于哈密顿函数
H
H
H，确定Metropolis准则，确定候选状态是否接受，如果
(
z
,
r
)
(\bm z, \bm r)
(z,r)是初始状态，
(
z
∗
,
r
∗
)
(\bm z^*,\bm r^*)
(z∗,r∗)是终止状态，那么接受概率为
如果蛙跳算法完美模拟哈密顿动力系统，那么
H
H
H值是不变的，但是由于数值误差，
H
H
H有时会减小，所以采用Metropolis准则消除这种偏差。这里
ϵ
\epsilon
ϵ要足够小。
- 这里的 ϵ \epsilon ϵ可以是前向，也可以是后向，可以以等概率随机选择
- 书上给的具有独立分量的多维高斯分布的例子， ϵ \epsilon ϵ的尺度受控于最小的标准差 σ min \sigma_{\min} σmin. 注意混合MC中蛙跳积分的目的是在相空间phase space中移动较大的距离到达新状态，这个新状态与初始状态相对独立，同时还有较高的接受率，为此必须多次蛙跳积分迭代，迭代次数是 σ max / σ min \sigma_{\max}/\sigma_{\min} σmax/σmin这一量级.
- 相比之下，在Metropolis算法中，随机游走得方式达到近似独立的状态则需要 ( σ max / σ min ) 2 (\sigma_{\max}/\sigma_{\min})^2 (σmax/σmin)2这一量级
11.6 Estimation the Partition Function
对于采样，很多方法绕过了
Z
E
Z_E
ZE。但它对于贝叶斯模型比较是有用的，因为表示了模型证据。
对于模型比较，需要计算配分函数之比，可以用如下重要性采样进行估计
其中
{
z
(
l
)
}
\{\bm z^{(l)}\}
{z(l)}从
p
G
(
z
)
p_G(\bm z)
pG(z)采样
- 如果
p
G
p_G
pG能很好地匹配
p
E
p_E
pE，那么这种方法会产生比较准确的结果。
如果两个分布差的太远，直接比不好比，可以在其中内插分布，得到 p 2 , … , p M − 1 p_2,\dots,p_{M-1} p2,…,pM−1
（这有点像花书上的退火重要采样）
一种确定中间分布的方法为构造如下渐变的能量函数
如果中间比值用MC算法得到，那么用一个单一的马尔可夫链要相对于每一个比值重新设置一个马尔可夫链的方式更高效。这种情况下，马尔可夫链初始设置为 p 1 p_1 p1，然后再某个合适的迭代次数后，移到序列中下一个概率分布。不过系统必须在每个阶段保持与均衡分布接近
另一种方法为从马尔可夫链中得到的样本来定义重要采样分布，样本集合为
z
(
1
)
,
…
,
z
(
L
)
\bm z^{(1)}, \dots, \bm z^{(L)}
z(1),…,z(L)
参考文献：
[1] Christopher M. Bishop. Pattern Recognition and Machine Learning. 2006
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
