PRML读书笔记——采样方法

本章中，我们希望解决的基本的问题涉及到关于⼀个概率分布p(z)寻找某个函数f(z)的期望。这⾥，z的元素可能是离散变量、连续变量或者⼆者的组合。因此，在连续变量的情形下，我们希望计算下⾯的期望

E[f]=∫f(z)p(z)dz

我们假设，使⽤解析的⽅法精确地求出这种期望是⼗分复杂的。因此，我们采用采样的方法抽取出样本，计算出对应的值

本章的采样方法主要包括：重要采样、拒绝采样、蒙特卡罗采样、吉布斯采样、切片采样和混合蒙特卡罗采样等

基本采样方法

这里，我们研究从一个给定的概率分布中生成随机样本的方法

标准概率分布

这考虑如何从简单的非均匀分布中生成随机数。可以从均匀分布开始，利用变换的方法，形成所需要计算的概率分布，并求出这个概率分布的不定积分的反函数。

假设z在区间(0,1)上均匀分布，我们使用某个函数f(.)对z进行变换，得到y=f(z)，则y上的概率分布为:

p(y)=p(z)∣∣∣dzdy∣∣∣

这里，y所具有的分布是我们希望的得到的分布形式。对上式积分，有：

z=h(y)≡∫y−∞p(yˆ)dyˆ

它就是p(y)的不定积分，因此y=h−1(z)。

这里存在两个难点：（1）变换f不容易构造（2）不定积分不一定容易求解

为了对非均匀分布采样，考虑另外的方法：拒绝采样和重要采样，但它们只能处理单变量概率分布的情况。

拒绝采样

拒绝采样的思想是，目前需要采样的概率分布形式复杂，我们引入相对简单的概率分布（提议分布，proposal distribution），这个概率分布能覆盖目标分布，然后直接在提议分布上采样，通过一定的准则选择拒绝还是接受当前的采样值。

假设我们希望从概率分布p(z)中采样，但是它不是简单的标准概率分布形式，从中采样困难。与大多数情况类似，假设对于任意给定的
z，我们能够很容易计算p(z)，即：

p(z)=1Zpp˜(z)

现在，我们引入简单的提议分布q(z)和常数k，使得：

∀z,kq(z)≥p˜(x)

则函数kq(z)被称为比较函数。整体情况下图说明。拒绝采样器的每个步骤涉及到⽣成两个随机数。⾸先，我们从概率分布q(z)中⽣成⼀个数z0。接下来，我们在区间[0,kq(z0)]上的均匀分布
中⽣成⼀个数u0。这对随机数在函数kq(z)的曲线下⽅是均匀分布。最后，如果u0>widetildep(z0)，那
么样本被拒绝，否则u0被保留。因此，如果它位于图11.4的灰⾊阴影部分，它就会被拒绝。这样，剩余的点对在曲线widetildep(z)下⽅是均匀分布的，因此对应的z值服从概率分布p(z)。

一个样本的接受率为p˜(z)kq(z)，因此整体样本的接受率为：

p(accepted)=∫{p˜(z)kq(z)}q(z)dz=1k∫p˜(z)dz

显然，k越大，整体接受率越低，因此k越低越好，但是k有需要满足限制kq(z)处处不小于p˜(z)

一般而言，q(z)的形式不啊后确定，可以直接基于概率分布p(z)构建提议分布的函数形式。首先，可以从ln p(z)的某些格点处开始计算，计算对应的切线，将各个切线连起来形成界限函数。然后从界限分布中采样，如果样本被接受，则他是所求概率分布的样本；反之，将它并入格点的机会中，计算新的切线，优化界限函数。随着格点数量的增加，界限函数对目标概率分布的近似效果逐渐变好。

---

缺点：因为目标概率分布的情况复杂，找到⼀个较好的提议分布和⽐较函数是⼀件相当困难的事情。此外，接受率随着维度的指数下降是拒绝采样的⼀个⼀般特征。虽然拒绝采样在⼀维或⼆维空间中是⼀个有⽤的⽅法，但是它不适⽤于⾼维空间

重要采样

这种采样方法主要用于估计概率分布的期望。与拒绝采样不同，重要采样不拒绝任何的采样结果，而是给提议分布上的采样结果赋予权重

同样假设直接从p(z)采样无法完成，但是给定z，p(z)很容易计算。为了计算期望，我们均匀地对z空间采样，然后计算期望：

E[f]=≃∑l=1Lp(z(l))f(z(l))

显然，这种采样是非常低小的，因为一般而言，目标概率分布都将它的大部分质量限制在z空间的一个很小的区域，也就是说只有很小部分的样本会对求合适产生贡献。

这里再次引入提议分布q(z)，期望可以表示为：

E(f)=∫f(z)p(z)dz=∫f(z)p(z)q(z)q(z)dz≃1L∑l=1Lp(z(l))q(z(l))f(z(l))

这里，rl=p(z(l))q(z(l))称为重要性权重，修正了由于从错误的概率分布中采样引⼊的偏差。注意，与拒绝采样不同，所有⽣成的样本都被保留。

考虑到归一化因子难以求解，因此，一般定义rl˜=p˜(z(l))q˜(z(l))

则：

E[f]≃∑l=1Lwlf(z(l))

wl=r˜l∑mr˜m=p˜(z(l))q(z(l))∑mp˜(z(m))q(z(m))

缺点：它具有产⽣任意错误的结果的可能性，并且这种错误⽆法检测。这是因为重要采样⽅法的成功严重依赖于采样分布q(z)与所求的概率分布p(z)的匹配程度。经常出现的情形是p(z)变化剧烈，并且⼤部分的质量集中于z空间的⼀个相对较⼩的区域中，此时重要性权重{rl}由⼏个具有较⼤值的权值控制，剩余的权值相对较⼩。因此，有效的样本集⼤⼩会⽐表⾯上的样本集⼤⼩L⼩得多。如果没有样本落在p(z)f(z)较⼤的区域中，那么问题会更加严重。此时，rl和rlf(z(l))的表⾯上的⽅差可能很⼩，即使期望的估计可能错得离谱。

采样-重要性-重采样

这个⽅法有两个阶段。在第⼀个阶段，L个样本z(1),...,z(L)从q(z)中抽取。然后在第⼆个阶段，用重要性采样的方法得到权重权值w1,...,wL。最后，L个样本的第⼆个集合从离散概率分布(z(1),...,z(L))中抽取，概率由权值(w1,...,wL)给定。

可以证明，⽣成的L个样本只是近似地服从p(z)，但是在极限L→∞的情况下，分布变为了正确的分布。

马尔科夫链蒙特卡罗（MCMC）

前面的基本采样方法在高维空间有很大的局限性。MCMC可以很好地应对样本空间维度的增长，MCMC的采样时针对时序状态的。

首先我们从提议分布中采样，记录下当前状态z(τ)，以及依赖于这个当前状态的提议分布q(z|z(τ))，从⽽样本序列z(1),z(2),...组成了⼀个马尔科夫链。如果我们有p(z)=p˜(z)Zp，那么假定对于任意的z值都可以计算p˜(z)，并且提议分布很简单，直接采样容易。在算法的每次迭代中，我们从提议分布中⽣成⼀个候选样本z∗，然后根据⼀个恰当的准则接受这个样本。

这里举一个例子，说明MCMC算法主要想解决的问题。考虑一个随机游走的例子。它有整数组成的状态空间z，概率为；

p(z(τ+1)=z(τ))=0.5

p(z(τ+1)=z(τ)+1)=0.25

p(z(τ+1)=z(τ)−1)=0.25

其中z(τ)表⽰在步骤τ的状态。如果初始状态是z(0)=0，那么根据对称性，在时刻τ的期望状态也是零，即E[z(τ)]=0，类似地很容易看到E[(z(τ))2]=τ/2 。因此，在τ步骤之后，随机游⾛所经过的平均距离正⽐于τ的平⽅根。这个平⽅根依赖关系是随机游⾛⾏为的⼀个典型性质，表明了随机游⾛在探索状态空间时是很低效的。设计马尔科夫链蒙特卡罗⽅法的⼀个中⼼⽬标就是避免随机游⾛⾏为。

马尔科夫链

这里主要分析马尔科夫链的基本性质——细节平衡性质（detailed balance）

首先，马尔科夫链具有齐次性，即下一时刻的状态只受当前状态的影响。

p(z(m+1)|z(1),...,z(m))=p(z(m+1)|z(m))

因此，边缘概率为：

p(z(m+1))=∑z(m)p(z(m+1)|z(m))p(z(m))

对于一个概率分布来说，如果马尔科夫链中的每⼀步都让这个概率分布保持不变，那么我们说这个概率分布关于这个马尔科夫链是不变的，或者静⽌的。对于⼀个转移概率为T(z′,z)的同质(对任意时刻，转移概率矩阵不变)的马尔科夫链来说，如果下式成立，则概率分布是不变的。

p∗(z)=∑z′T(z′,z)p∗(z′)

细节平衡性质：是确保所求的概率分布p(z)不变的⼀个充分（⾮必要）条件，如下式。

p∗(z)T(z,z′)=p∗(z′)T(z′,z)

满足细节平衡就能收敛到平稳分布，如下：

∑z′p∗(z′)T(z′,z)=∑z′p∗(z)T(z,z′)=p∗(z)∑z′p(z|z′)=p∗(z)

这里就解释细节平衡得必要性，如果构造的马尔科夫链是的所求的概率分布是不变的，我们就能控制它收敛到指定的分布p∗(z)。

Metropolis 算法

我们假定提议分布是对称的，即q(zA|zB)=q(zB|zA)对于所有的zA和zB成⽴。这样，候选的样本被接受的概率为:

A(z∗,z(τ))=min(1,p˜(z∗)p˜(z(τ)))

可以这样实现：在单位区间(0, 1)上的均匀分布中随机选择⼀个数u，然后如果A(z∗,z(τ))>u就接受这个样本。注意，如果从z(τ) 到z∗引起了p(z)的值的增⼤，那么这个候选样本会被保留。

如果候选样本被接受，那么z(τ+1)=z∗，否则候选样本点z∗被丢弃，z(τ+1)被设置为z(τ)，
然后从概率分布q(z|z(τ+1))中再次抽取⼀个候选样本。当⼀个候选点被拒绝时，前⼀个样本点会被包含到最终的样本的列表中（不会被直接丢弃），从⽽产⽣了样本点的多个副本。

Metropolis-Hastings 算法

它是基于metropolis的改进，这里提议分布不再是对称的，而且接受概率的形式发生改变。在算法的步骤τ中，当前状态为z(τ)，我们从概率分布qk(z|z(τ))中抽取⼀个样本z∗，则接受率为：

A(z∗,z(τ))=min(1,p˜(z∗)qk(z(τ)|z∗)p˜(z(τ))qk(z∗|z(τ)))

这样设计的理由是，形成的马尔科夫链是满足细节平衡的，也即是不变的马尔科夫链。即使随机游走，也能确保收敛到平稳分布。

吉布斯采样

吉布斯采样是马尔科夫链蒙特卡罗算法的具体情形

考虑我们想采样的概率分布p(z)=p(z1,...,zM)，并且假设我们已经选择了马尔科夫链的某个初始的状态。吉布斯采样的每个步骤涉及到将⼀个变量的值替换为以剩余变量的值为条件，从这个概率分布中抽取的那个变量的值

一般步骤如下：

（1）初始化{zi:i=1,...,M}

（2）对于τ=1,...,T：

• 采样 zτ+11∼p(z1|zτ2,zτ3,...zτM)

• 采样 zτ+12∼p(z2|zτ+11,zτ3,...zτM)

… …

• 采样 zτ+1j∼p(zj|zτ+11,zτ+12,...zτ+1j−1,zτj+1,...,zτM)

… …

• 采样 zτ+1M∼p(zM|zτ+11,zτ+12,...zτ+1M−1)

切片采样

因为Metropolis算法对步长敏感。步长小，会产生随机游走，算法很慢；步长大，拒绝率较高，算法低效。切片采样替工了可以自动调节步长来匹配分布特征的方法。

这种采样方法使用额外的变量u对z进行增广，然后从联合的(z,u)空间中采样。比如，概率分布如下：

pˆ(z,u)=⎧⎩⎨⎪⎪1Zp00≤u≤p˜(z)others

则z的边缘概率分布为：

∫pˆ(z,u)du=∫p˜(z)01Zpdu=p˜(z)Zp=p(z)

我们可以通过从pˆ(z,u)中采样，然后忽略u值的⽅式得到p(z)的样本。通过交替地对z和u进⾏采样即可完成这⼀点。给定z的值，我们可以计算p˜(z)的值，然后在0≤u≤p˜(z)上均匀地对u进⾏采样，这很容易。然后，我们固定u，在由{z:p˜(z)>u}定义的分布的“切⽚”上，对z进⾏均匀地采样。

显然，得到u后，直接从穿过概率分布的切⽚中采样很困难。近似方法如下：

从⼀个包含z(τ)的具有某个宽度w的区域开始，然后测试每个端点，看它们是否位于切⽚内部。如果有端点没在切⽚内部，那么区域在增加w值的⽅向上进⾏扩展，知道端点位于区域外。然后，z′的⼀个样本被从这个区域中均匀抽取。如果它位于切⽚内，那么它就构成了z(τ+1)。如果它位于切⽚外，那么区域收缩，使得z′组成⼀个端点，并且区域仍然包含z(τ)。然后，另⼀个样本点从这个缩⼩的区域中均匀抽取，以此类推，直到找到位于切⽚内部的⼀个z值。

混合蒙特卡罗采样

Metroplolis算法的主要局限是受步长影响，具有随机游走的行为。混合蒙特卡罗算法将哈密顿动态系统与Metropolis算法相结合，消除与离散化过程关联的任何偏差。

算法使⽤了⼀个马尔科夫链，它由对动量变量r的随机更新以及使⽤蛙跳算法对哈密顿动态系统的更新交替组成。在每次应⽤蛙跳算法之后，基于哈密顿函数H的值，确定Metropolis准则，确定⽣成的候选状态被接受或者拒绝。因此，如果(z,r)是初始状态，(z∗,r∗)是蛙跳积分后的状态，那么候选状态被接受的概率为

min(1,exp{H(x,r)−H(z∗,r∗)})

这里的H就是哈密顿动态系统的总能量，为势能和动能之和。感觉这部分没看懂，所以具体推导和分析就省略了