1. 马尔可夫链

1.1. 马尔可夫性质 & 转移矩阵

马尔可夫链是一个随机过程 ${x_i}$ ，它满足以下条件： $$P(x_i|x_0,\ldots,x_{i-1})=p(x_i|x_{i-1})$$ 我们可以通俗地理解为，系统在未来时刻 $i$ 的状态只取决于它当前时刻 $i−1$ 的状态，而与它过去的所有状态 $(x_0,\ldots,x_{i-2})$ 都无关 。

我们定义随机过程 ${x_i}$ 所有可能取值的集合为状态空间。状态空间可以是有限的，也可以是无限的。对于状态有限的离散马尔可夫链，我们可以用一个矩阵 $P$ 来描述状态之间转移的规则 。

矩阵中的每一个元素 $p_{ij}​$ 代表了从状态 $i$ 转移到状态 $j$ 的概率 ： $$p_{ij}=p(x_n=j|x_{n-1}=i)$$ 这个矩阵有一个很重要的性质，即： $$\sum_jp_{ij}=1$$ 因为从任何一个状态 $i$ 出发，它下一步必须转移到状态空间中的某一个状态（包括可能停在原地），所以所有可能转移的概率之和必然等于1。一个3x3的转移矩阵例子如下： $$P=\begin{pmatrix}\frac{1}{2}&\frac{1}{2}&0\\frac{1}{2}&0&\frac{1}{2}\0&\frac{1}{2}&\frac{1}{2}\end{pmatrix}$$

我们可以这样解读这个矩阵：

• 第一行: 如果当前在状态1，那么下一步有1/2的概率留在状态1，有1/2的概率转移到状态2，有0的概率直接转移到状态3。
• 第二行: 如果当前在状态2，那么下一步有1/2的概率转移到状态1，有0的概率留在状态2，有1/2的概率转移到状态3。
• 第三行: 如果当前在状态3，那么下一步有1/2的概率转移到状态2，有1/2的概率留在状态3，有0的概率直接转移到状态1。

1.2. 正则矩阵

如果存在一个正整数 $k$，使得矩阵 $P^k$ 的所有元素都严格为正（即大于0），那么我们称转移矩阵 $P$ 是正则的 (regular) 。“正则性”意味着，从任何一个状态出发，经过固定的 $k$ 步之后，都有可能到达任何其他状态。

$P^k$ 的元素 $(P^k)_{ij}$ 代表的就是从状态 $i$ 经过 $k$ 步恰好到达状态 $j$ 的概率。所有元素大于0，就保证了任意两个状态之间在 $k$ 步后都是连通的。

一个马尔可夫链如果其转移矩阵是正则的，那么它就具有一个非常好的性质：无论从哪个状态开始，经过足够长的时间后，它都会收敛到一个唯一的平稳分布。这是MCMC方法能够工作的理论基石之一。

计算一下上面例子中矩阵的平方 $P^2$： $$P^2=P\times P=\begin{pmatrix}\frac{1}{2}&\frac{1}{2}&0\\frac{1}{2}&0&\frac{1}{2}\0&\frac{1}{2}&\frac{1}{2}\end{pmatrix}\begin{pmatrix}\frac{1}{2}&\frac{1}{2}&0\\frac{1}{2}&0&\frac{1}{2}\0&\frac{1}{2}&\frac{1}{2}\end{pmatrix}=\begin{pmatrix}\frac{1}{2}&\frac{1}{4}&\frac{1}{4}\\frac{1}{4}&\frac{1}{2}&\frac{1}{4}\\frac{1}{4}&\frac{1}{4}&\frac{1}{2}\end{pmatrix}$$ $P^2$ 的所有元素都大于0。因此，这个例子中的转移矩阵 $P$ 是一个正则矩阵（取 $k=2$ 即可）。这意味着，无论你从状态1、2还是3开始，走2步之后，你都有可能出现在任何一个状态。

从状态 $i$ 出发，经过两步之后到达状态 $j$ 的概率被记为 $p_{ij}^{(2)}$ 。它的计算公式是： $$p_{ij}^{(2)}=p(x_{n+1}=j|x_{n-1}=i)=\sum_kp_{ik}p_{kj}=(P^2){ij}$$ 从状态 $i$ 要走两步到状态 $j$，它必须在第一步时先到达一个中间状态，我们称之为 $k$。从 $i$ 到 $k$ 的概率是 $p{ik}$，从 $k$ 到 $j$ 的概率是 $p_{kj}​$，所以，通过特定中间状态 $k$ 完成这个两步转移的概率就是 $p_{ik}\times p_{kj}$ 。这个中间状态 $k$ 可以是状态空间里的任何一个状态。因此，我们需要把所有可能的中间路径的概率全部加起来，就得到了从 $i$ 两步到 $j$ 的总概率，因此要求和。这个求和公式 $\sum_kp_{ik}p_{kj}$ 正好是矩阵乘法的定义。

思路延伸下去，从初始状态（时刻0）为 $i$ 出发，经过 $n$ 步后到达状态 $j$ 的概率，其对应的矩阵就是 $P$ 的 $n$ 次方。 $$p(x_n=j|x_0=i)=p_{ij}^{(n)}$$ 而包含所有这些 $n$ 步转移概率的矩阵就是： $$P^n=P\times P\times\cdots\times P\quad(n\mathrm{~times})$$ 这个性质是马尔可夫链理论中非常有用的一个工具。它意味着只要我们知道了一步转移矩阵 $P$，我们就可以通过计算矩阵的幂来预测系统在未来任意时刻的状态概率分布。

2. 平稳链(Stationary Chain)

2.1. 平稳分布 (stationary distribution)

平稳分布 $π$ 的数学定义：如果一个分布 $π$ 满足以下方程，那么它就是该马尔可夫链的平稳分布。 矩阵形式: $$\pi^\top=\pi^\top P$$ 元素形式: $$\pi_j=\sum_i\pi_ip_{ij}$$ 如果当前时刻系统的状态分布是 $π$ (一个描述每个状态被占据的概率的行向量)，那么经过一次转移（即乘以转移矩阵 $P$ ）之后，系统在下一时刻的状态分布依然是 $π$。分布不再随时间发生变化，达到了一个“平稳”或者“均衡”的状态。

从线性代数的角度看，这个公式 $\pi^\top=\pi^\top P$ 表明，平稳分布 $\pi^\top$ 是转移矩阵 $P$ 的一个左特征向量，其对应的特征值为1。

假设我们的初始状态 $x_0$​ 就是从平稳分布 $π$ 中抽取的，也就是说 $p(x_0=i)=\pi_i$ ，那么在下一时刻 $x_1​$ 系统处于状态 $j$ 的概率为： $$p(x_1=j)=\sum_ip(x_1=j|x_0=i)p(x_0=i)=\pi_j$$ 这个推导表明，如果系统以平稳分布开始，那么在下一步之后，系统的分布仍然是这个平稳分布。通过数学归纳法，可以证明这个性质对于任意步数 $n$ 都成立。也就是说，一旦马尔可夫链的分布进入了平稳分布 $π$，它将永远保持在这个分布上。

2.2. 平稳分布的存在性

结论：如果一个马尔可夫链的转移矩阵 $P$ 是正则的，那么它必然拥有一个极限概率分布 $(\pi_1,\ldots,\pi_N)$，这个极限分布 $π$ 具备两个性质：

• 它是一个合法的概率分布，即所有状态的概率之和为1: $\sum_j\pi_j=1$。
• $n$ 步转移概率的极限等于这个分布的对应分量： $$\lim_{n\to\infty}p_{ij}^{(n)}=\pi_j$$ 这告诉我们，极限分布与初始状态无关。即无论从哪里开始，长期来看，马尔可夫链都会收敛到 $π$。在很长一段时间后，我们在状态 $j$ 找到这个链的概率约等于 $π_j​$，这个概率不依赖于它的起始点。

这个“遗忘初始状态”的特性，是马尔可夫链蒙特卡洛（MCMC）方法能够成功应用的核心原因。在解决实际问题时（比如在贝叶斯统计中），我们想要从一个非常复杂的目标分布（比如后验分布）中抽样。MCMC的策略就是巧妙地设计一个马尔可夫链，使得这个链的平稳分布正好就是我们想要的那个目标分布。

我们不需要知道如何直接从这个复杂的目标分布开始抽样。我们可以从一个任意的、随便选择的初始点开始运行我们设计的马尔可夫链。只要我们让链运行足够长的时间（这个过程通常被称为“预烧/老化 (burn-in)”），它就会逐渐“忘记”它那个随便选择的起点，并最终收敛到其唯一的平稳分布，也就是我们的目标分布。

2.3. 平稳分布的唯一性

结论：对于一个在 $N$ 个状态上定义的正则 转移矩阵 $P$，其极限（也就是平稳）分布 $π$ 是下面这个方程组的唯一的非负解：

1. $\pi_j=\sum_{k=1}^N\pi_kp_{kj}$ （矩阵形式即 $\pi^\top=\pi^\top P$）
2. $\sum_{k=1}^N\pi_k=1$

这个性质把寻找平稳分布的问题，转化成了一个具体的代数问题：求解一个有 $N+1$ 方程的线性方程组（$N$ 个来自 $\pi^\top=\pi^\top P$，1个来自概率求和为1）。只要解出这个方程组，我们就能精确地知道链的长期行为了。

2.4. 平稳分布的物理解释：时间占比

如果 $n$ 步转移概率 $p_{ij}^{(n)}$ 的极限存在且等于 $π_j​$，那么这些概率在时间上的平均值也会收敛到同一个极限。 $$\text{如果 }\lim_{n\to\infty}p_{ij}^{(n)}=\pi_j,\quad\text{那么 }\lim_{m\to\infty}\frac{1}{m}\sum_{k=1}^mp_{ij}^{(k)}=\pi_j$$ 这个“平均概率”的内在含义如下： $$\frac{1}{m}\sum_{k=1}^mp_{ij}^{(k)}=\frac{1}{m}E[\sum_{k=1}^mI(x_k=j|x_0=i)]$$

• $I(x_k=j)$ 是一个指示函数：如果在第 $k$ 步时，链处于状态 $j$，那么它的值为1，否则为0。
• $\sum_{k=1}^mI(x_k=j)$ 就是在总共 m 步的行程中，链访问状态 j 的总次数。
• $\frac{1}{m}\sum_{k=1}^mI(x_k=j)$ 就是链在 m 步中，停留在状态 j 的时间比例（或者说频率）。
• $E[…]$ 表示对这个比例取期望值，即平均时间比例。

将所有这些联系起来，我们得到了一个非常直观且重要的结论：$π_j​$ 就等于马尔可夫链在长期运行中，停留在状态 j 的平均时间比例。

平稳分布 $π_j​$ 具有双重含义：

1. 概率极限: 无论从哪个状态出发，经过无穷多步后，系统处于状态 $j$ 的概率就是 $π_j​$。
2. 时间占比: 在一个长期运行的链中，系统停留在状态 $j$ 的时间约占总时间的 $π_j​$ 比例。

例如，一个天气的马尔可夫模型，如果我们计算出其平稳分布中“雨天”对应的分量 $π_{雨天}​=0.1$，这就意味着：

1. 从长期来看，未来任意一天是雨天的概率是10%。
2. 在过去很长一段时间里（比如1000天），大约有10%的日子（100天）是雨天。

2.5. 遍历性 (Ergodic)

一个正则 的马尔可夫链具有非常好的性质，比如存在唯一的平稳分布。而保证马尔可夫链具有良好长期行为的三个核心属性，这三个属性共同构成了遍历性 (Ergodicity)。

1. 不可约性 (Irreducible)

定义: 对于状态空间中的任意两个状态 $(i,j)$，都存在一个步数 $n>0$，使得从状态 $i$ 经过 $n$ 步转移到状态 $j$ 的概率大于零，即 $p_{ij}^{(n)}>0$ 。

不可约性保证了整个状态空间是完全连通的。可以把状态想象成城市，把转移概率想象成道路。不可约就意味着，从任何一个城市出发，总能找到一条路到达任何其他城市。

在MCMC中，如果链不是不可约的，就意味着有些状态可能永远也访问不到，那我们采样的样本就无法代表整个目标分布，得到的结果将是有偏的。

2. 非周期性 (Aperiodic)

定义: 链不会以固定的、规律性的时间间隔返回到某个状态子集。

这个性质是为了排除一种特殊情况：链被困在一个确定性的循环里。在MCMC中，我们希望链能够在状态空间中自由地探索，而不是被一个固定的“节拍”束缚住。如果一个链是周期的，它的收敛可能会出现震荡，而不是平稳地趋近于平稳分布。通常，只要链中有一个状态存在自环（即 $p_{ii}​>0$），就可以破坏周期性。

3. 常返性 (Recurrent)

对于任意给定的状态 $i$，如果链从状态 $i$ 出发，它最终能返回到状态 $i$ 的概率为1。常返性保证了链不会“一去不复返”。无论它游走到多远的地方，最终总会回到它曾经访问过的任何一个状态。

• 正常返 (Positive Recurrent): 不仅保证能返回，而且返回的平均时间是有限的。这意味着链不会离家太久。对于有限状态的马尔可夫链，只要是不可约的，就一定是正常返的。
• 零常返 (Null Recurrent): 保证能返回，但返回的平均时间是无限的。这通常发生在某些无限状态空间的链中，链可能会进行极长时间的“远足”才回家。

只有正常返的链才能保证平稳分布的存在。链需要频繁地返回各个状态，才能让长期的平均行为稳定下来。

当一个马尔可夫链同时满足不可约、非周期、且正常返这三个条件时，我们就称这个链是遍历的 (ergodic)。一个遍历的马尔可夫链，是我们在MCMC中最理想的模型。它保证了链不仅存在唯一的平稳分布，而且无论从哪里开始，链的长期行为都会收敛到这个平稳分布。

对于遍历的马尔可夫链，有一个强大的遍历定理。它表明，当模拟时间足够长时，对某个函数在链的轨迹上计算的时间平均，会收敛到该函数在平稳分布下的空间平均（即期望）。

这正是我们用MCMC来计算复杂积分的理论依据：我们通过模拟一个很长的链（计算时间平均），来近似我们无法直接计算的期望值（空间平均）。

3. 正则转移概率矩阵

定理：对于一个有限状态的马尔可夫链，它是非周期性且不可约的当且仅当它的转移概率矩阵是正则的。基于这个等价关系，可以推论：极限分布 $π$ 是存在的。

一个链要具备良好的长期收敛性质，最好是“遍历的”（即不可约、非周期、正常返）。对于有限状态的链，只要满足不可约和非周期这两个条件，它就是遍历的。

但直接从定义去判断“不可约”和“非周期”有时会比较抽象。而判断一个矩阵是否“正则”则是一个非常具体的操作：我们只需要计算它的幂 $P^2,P^3,\ldots$ ，看看是否存在某一个 $P^k$ 它的所有元素都大于零。

4. 可逆性 (reversibility)

对于一个非周期、不可约 的马尔可夫链：

• 它的平稳分布 $π$ 是唯一的。
• $n$ 步转移矩阵 $P_n$ 会收敛到一个特殊的矩阵。 $$\lim_{n\to\infty}P^n=\pi^\top$$ 它意味着，当 $n$ 趋于无穷大时，转移矩阵 $P_n$ 会变成一个每一行都完全相同的矩阵，而它的每一行都是平稳分布向量 $\pi^\top$。

当步数 $n$ 足够大时，我们从链中得到的当前状态 $x^{(n)}$，就可以被看作是直接从平稳分布 $π$ 中抽取的一个样本。我们设计一个以目标分布为平稳分布的链，让它运行足够长的时间（burn-in），然后收集它产生的状态，这些状态就是我们想要的样本。

果一个马尔可夫链的转移概率和平稳分布满足以下方程，我们就称这个链是可逆的 (reversible)。 $$\pi_ip_{ij}=\pi_jp_{ji}$$ 这个方程有一个更常用的名字，叫做细致平稳条件。这个公式的含义是，在平稳状态下，从状态 $i$ 离开并到达状态 $j$ 的“流量”，与从状态 $j$ 离开并到达状态 $i$ 的“流量正好相等。

5. Metropolis-Hastings (M-H) 算法

5.1. 思想

在统计和机器学习中，我们经常需要计算某个函数 $h(x)$ 在概率分布 $f(x)$ 下的期望值，也就是积分 $\int h(x)f(x)dx$ 。比如，在贝叶斯统计中，$f(x)$ 可能是一个非常复杂的后验分布，而我们想计算这个分布的均值、方差等统计量。

解决这个积分的传统方法是蒙特卡洛法：从 $f(x)$ 中抽取大量独立样本 $X_1,...,X_n$ 然后用均值 $\frac{1}{n}\sum h(X_i)$ 来近似积分值。但最大的障碍是：$f(x)$ 的形式通常非常复杂，我们根本没有办法直接从中抽取独立的样本。

我们可以通过构建一个遍历的 (ergodic) 马尔可夫链，来获得一系列样本 $X_1,...,X_n$，这些样本可以用来近似 $f(x)$，从而避免了直接从 $f(x)$ 抽样的难题。

因为马尔可夫链具有无记忆性 $P(x_i|x_0,\ldots,x_{i-1})=p(x_i|x_{i-1})$ ，我们可以从一个任意的初始点 $x^{(0)}$ 开始，使用一个特定的转移核（也就是M-H算法本身）来生成一个马尔可夫链，这个转移核被精心设计，以确保其平稳分布恰好是 $f(x)$。

只要时间足够长，在各个状态的停留比例最终都会趋向于那个由$f(x)$决定的唯一平稳分布。

5.2. MCMC

• 任何一个能够生成一个遍历的 (ergodic) 马尔可夫链 ($X^{(t)}$)，并且该链的平稳分布恰好是我们想要模拟的目标分布 $f(x)$ 的方法，都称之为MCMC方法。实现MCMC的两个关键概念呢如下：

目标密度 $f(x)$:

• 这就是我们最终的目标。它通常形式复杂、维度很高（例如贝叶斯后验分布），导致我们无法直接从中采样。

提议/候选密度 $q(x,y)$ 或 $q(y∣x)$:

• $q(y|x)$ 的含义是：给定我们当前的位置在 x，我们提议下一步移动到 y 的概率分布是什么。这个提议分布是我们自己选择的，关键是要容易从中采样。比如，一个常见的选择是以当前点 $x$ 为中心的正态分布。我们每走一步，都从这个正态分布中抽个样，得到一个“提议”的新位置。

5.3. 具体执行

假设我们当前处在链的第 $n-1$ 步，状态为 $x_{n−1}​=x$。要生成下一步的状态 $x_n$​，我们执行以下操作：

1. 提议 (Proposal)

• 从我们预先选择的提议分布 $q(y|x)$中抽取一个候选样本 $y$。
• 这一步就像我们在“藏宝图”上，根据一个简单的移动规则（比如，以当前位置为中心，随机看一个方向和距离），提出了一个想要去的新地点 $y$。

2. 决策 (Decision)

• 我们计算一个接受概率 $α(x,y)$，并以这个概率来决定是否接受候选样本 $y$。
• 接受概率的公式: $$\alpha(x,y)=\min\left(\frac{f(y)q(y,x)}{f(x)q(x,y)},1\right)$$ 这个公式是M-H算法的精髓，我们可以把它拆成两部分来理解：
• 目标比率$\frac{f(y)}{f(x)}$: 这一部分比较的是新提议的地点 $y$ 与当前地点 $x$ 在“藏宝图”上的“价值”。如果 $f(y)>f(x)$，说明新地点比当前地点更好（概率密度更高），这个比率就大于1，使得我们倾向于接受这个移动。如果 $f(y)<f(x)$，则新地点更“差”，比率小于1，我们就以一定概率接受这个“更差”的移动。允许向“更差”的地方移动是算法能够跳出局部最优，探索整个分布的关键。
• 提议比率 (修正项) $\frac{q(y,x)}{q(x,y)}$​: 这是Hastings引入的修正项，它用来修正提议分布的不对称性。如果从 $x$ 跳到 $y$ 的概率 ($q(y∣x)$) 比从 $y$ 跳回 $x$ 的概率 ($q(x∣y)$) 更大，那么这个修正项就会变小，从而降低接受概率。这可以防止我们的链偏向于那些“易进难出”的区域，保证了探索的公平性。这个修正项是确保算法满足细致平稳条件 的关键。
• $min(..., 1)$ 的作用是确保概率值不会超过1。

我们通常是生成一个在 $(0, 1)$ 之间均匀分布的随机数 $u$，如果 $u≤α(x,y)$，我们就接受提议。

3. 更新 (Update)

• 规则: 如果提议 $y$ 被接受，那么链的下一个状态就是 $x_n​=y$；否则，如果提议被拒绝，链就停在原地不动，下一个状态仍然是 $x_n​=x$。
  • 拒绝样本的作用: 与拒绝采样不同，MCMC中的“拒绝”不是把样本丢掉。当一个提议被拒绝时，我们会将当前样本再次记录到我们的样本序列中。这非常重要，因为它使得链在高概率区域停留的时间更长（因为从高概率区向低概率区的移动更容易被拒绝），从而自动地让采集到的样本符合目标分布 $f(x)$ 的疏密程度。

M-H算法在计算接受概率的比率时，目标分布 $f(x)$ 的归一化常数（通常很难计算）被消掉了。我们只需要知道 $f(x)$ 的函数形式即可。它可以与任何方便采样的提议分布 $q(y∣x)$ 结合使用，算法的修正项会自动校正其中的偏差。

5.4. 数学化描述

从第 i 步到第 i+1 步的完整转移过程。

• 初始化: 算法从一个任意选择的初始值 $x^{(0)}$ 开始。
• 在第 i 步: 假设我们当前的状态是 $x^{(i)}$，我们执行以下两步：
  1. 提议 (Proposal): 从提议分布 $q(y∣x^{(i)})$ 中生成一个候选的新状态 $y$。
  2. 更新 (Update): 在第 $i+1$ 步，我们根据一个概率规则来设置新状态 $x^{(i+1)}$。这个规则可以写成一个简洁的分段函数： $$x^{(i+1)}=\begin{cases}y,&\text{以概率 }\alpha(x^{(i)},y)\text{ 发生}\x^{(i)},&\text{以概率 }1-\alpha(x^{(i)},y)\text{ 发生}&&\end{cases}$$ 这里的 $α(x^{(i)},y)$ 就是我们之前讨论过的接受概率。计算公式如下： $$\alpha(x,y)=\min\left{\frac{f(y)q(x|y)}{f(x)q(y|x)},1\right}$$
• 这里的 $q(y∣x)$ 表示“给定当前状态是x，提议下一步去y的概率”，而 $q(x∣y)$ 则是“给定状态是y，提议下一步去x的概率”。
• 这种写法比之前的 $q(x,y)$ 更能清晰地体现出转移的方向性，也让我们能更清楚地理解修正项 $q(y∣x)/q(x∣y)​$ 的含义：它是反向转移提议概率与正向转移提议概率的比值。

5.5. 算法特点

5.5.1. 接受“更好”的提议

- **倾向于“向上爬”**: 如果一个提议的移动使得接受概率公式中的比率 $\frac{f(y)q(x|y)}{f(x)q(y|x)}$​ 大于1，那么这个提议**总是会被接受**（接受概率为1）。这通常发生在当我们从一个低概率区域移动到一个高概率区域时。

- **也可能“向下走”**: 即使一个提议的移动使得上述比率小于1（比如移动到概率更低的区域），算法**仍然有可能接受它**。

• 登山者比喻:
  • 我们可以把M-H算法想象成一个在雾中探索山脉（目标分布）的“醉醺醺的登山者”。
  • 他倾向于往高处走（“向上爬”），因为这样更有可能接近山顶（高概率区域）。
  • 但为了不被困在某个小山头（局部最优解）上，他又会以一定的概率选择往下走一步（“向下走”）。正是这种偶尔“向下走”的能力，保证了他能够探索整个山脉，而不是只停留在第一个发现的山峰上。这与随机优化 的思想非常相似。

5.5.2. 无需归一化常数

• M-H算法一个巨大的实践优势：它只依赖于概率密度的比率，例如 $f(y)/f(x)$，因此与归一化常数无关。
  • 在许多实际问题中（尤其是贝叶斯推断），我们的目标分布 $f(x)$ 通常只能写成 $f(x)=\frac{1}{Z}\tilde{f}(x)$ 的形式。其中，$\tilde{f}(x)$ 的函数形式我们是知道的（比如 正比于 似然 × 先验），但归一化常数 $Z$ 是一个非常复杂、难以计算甚至无法计算的积分。
  • M-H算法的优势在于，当我们计算接受概率中的核心比率时，$Z$ 会被完美地消掉：$\frac{f(y)}{f(x)}=\frac{\tilde{f}(y)/Z}{\tilde{f}(x)/Z}=\frac{\tilde{f}(y)}{\tilde{f}(x)}$
  • 这意味着，我们不需要知道 $Z$ 到底是多少，就可以运行整个算法。这是M-H算法能够被广泛应用的最主要原因之一。

5.5.3. MCMC样本的特性

• M-H算法生成的样本与我们通常所说的独立同分布 (i.i.d.) 样本有本质区别。
  • 1. 样本是相关的: MCMC生成的是一个马尔可夫链，链中的每一个样本 $x^{(t)}$ 都是从前一个样本 $x^{(t−1)}$ 转移而来的，所以它们之间存在自相关性。
  • 2. 样本会重复: 当一个提议被拒绝时，算法并不会丢弃这个样本，而是将当前状态再次加入到样本序列中。因此，在MCMC的输出序列中看到连续出现相同的数值是非常正常的。
  • 重复的意义: 这些重复值并不是无用的。一个状态被重复的次数，恰恰反映了它在目标分布 $f(x)$ 中所占的“权重”。高概率区域的点更难“跳出去”（因为向外跳的提议更容易被拒绝），因此会被重复更多次，从而使得最终样本的直方图能够准确地还原 $f(x)$ 的形状。

5.6. 随机游走Metropolis算法

Random Walk Metropolis Algorithm是Metropolis-Hastings算法中一个非常重要且常见的特例。它的核心思想是通过选择一种特殊的、对称的提议分布，来大大简化计算。

这种算法的关键在于算法的提议分布 $q(x,y)$。它的形式是 $q(x,y)=g(y−x)$，其中 $g$ 是一个球对称 的函数。“对称”的直观含义是，从 $x$ 跳到 $y$ 的提议概率，和从 $y$ 跳回到 $x$ 的提议概率是完全相等的。也就是说，$q(y∣x)=q(x∣y)$。

实现这种对称提议的一个非常自然的方式是：$y=x+\epsilon$ 。这里的 $ϵ$ 是从一个关于原点对称的分布 $g$ 中抽取的随机扰动（或称为“步长”）。这个过程就像是在当前位置 $x$ 的基础上，进行一次“随机游走”。这就是这个算法名字的由来。

关于原点对称的分布 $g$ 有很多常见的选择，比如：

• 正态分布 (高斯分布): $ϵ∼N(0,σ^2)$。这是最常用的选择。
• 均匀分布: $ϵ∼Unif[−a,a]$。

从 $x$ 提议 $y$ 是基于扰动 $ϵ=y−x$ 发生的，其概率密度为 $g(ϵ)$。而从 $y$ 提议 $x$ 则是基于扰动 $x−y=−ϵ$ 发生的，其概率密度为 $g(−ϵ)$。因为 $g$ 是对称的，所以 $g(ϵ)=g(−ϵ)$，这就保证了 $q(y∣x)=q(x∣y)$。

因为设计了对称的提议分布，所以一般M-H算法接受概率公式中的修正项（提议比率）就等于1了。 $$\alpha(x,y)=\min\left(\frac{f(y)q(y,x)}{f(x)q(x,y)},1\right)=\min\left(\frac{f(y)}{f(x)},1\right)$$ 所以步骤如下：

从当前状态 $x_{n−1}​=x$ 开始：

1. 提议 (Proposal):
   • 从一个对称分布 $g$ 中（如 $N(0,σ^2)$）抽取一个随机扰动 $ϵ$。
   • 生成候选点 $y=x+ϵ$。
2. 决策 (Decision):
   • 生成一个从0到1的均匀分布随机数 $u∼Unif[0,1]$。
3. 更新 (Update):
   • 判断是否 $u\leq\min(\frac{f(y)}{f(x)},1)$。
   • 如果成立，则接受提议，令下一个状态 $x_n​=y$。
   • 如果不成立，则拒绝提议，令下一个状态保持不变 $x_n​=x$。

这个简化后的算法直观且易于实现，是MCMC方法中最基础和最常用的构建模块之一。我们来看一个例子：

我们的目标是生成一个服从标准正态分布的随机变量。其概率密度函数 (PDF) 通常用 $ϕ(⋅)$ 表示。 $$f(x)=\phi(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$$ 我们选择均匀分布作为我们的提议分布 $g$。具体来说，我们每次的随机“步长” $ϵ$ 是从一个对称的均匀分布 $Unif[−δ,δ]$ 中抽取的。这里的 $δ$ 是一个需要我们自己设定的参数。

• 提议: 从均匀分布 $Unif[−δ,δ]$ 中随机抽取一个步长 ϵ。
• 生成候选点: 候选点 $y=x+ϵ$，其中 $x$ 是当前位置。
• 决策: 生成一个0到1之间的随机数 $u$。
• 更新: 计算接受概率，并判断是否接受。由于提议分布是对称的，我们使用简化的Metropolis接受准则：
• 如果 $u\leq\min\left(\frac{\phi(y)}{\phi(x)},1\right)$，则接受提议，移动到新位置。
• 补充 (具体计算): 在实际编程中，这个比率计算起来很简单： $$\frac{\phi(y)}{\phi(x)}=\frac{\frac{1}{\sqrt{2\pi}}e^{-y^2/2}}{\frac{1}{\sqrt{2\pi}}e^{-x^2/2}}=e^{-(y^2-x^2)/2}$$ 参数 $δ$ 控制了我们“随机游走”的步长范围，它的取值对算法的效率有决定性的影响。
• 如果 $δ$ 太小 (例如 $δ=0.1$):
  • 现象: 每次提议的步子都非常小，候选点 $y$ 总是离当前点 $x$ 很近。
  • 结果: $ϕ(y)$ 和 $ϕ(x)$ 的值会非常接近，导致比率 $\frac{\phi(y)}{\phi(x)}$ 几乎总是约等于1。因此，接受率会非常高。
  • 问题: 链会像一个步履蹒跚的老人，移动得极其缓慢，需要非常非常多步才能探索整个分布。
• 如果 $δ$ 太大 (例如 $δ=10$):
  • 现象: 每次提议的步子都非常大，相当于在地图上进行“乾坤大挪移”。
  • 结果: 如果当前位置 $x$ 在概率高的中心区域（比如0附近），那么一个大步很有可能直接跳到概率极低的尾部区域。这将导致比率 $\frac{\phi(y)}{\phi(x)}$​ 非常小。因此，接受率会非常低。
  • 问题: 链会在绝大多数时间里都“原地踏步”（因为提议总被拒绝），偶尔才成功跳一步。探索效率同样低下。

这里的问题是：高自相关意味着差的混合。

自相关: 指的是MCMC样本序列中，一个样本 $x^{(t)}$ 与其后续样本 $x^{(t+k)}$ 之间的相关性。由于链的构造方式，样本之间必然是相关的。

混合: 指的是马尔可夫链探索整个目标分布的速度和效率。一个“混合良好”的链，能够快速地在分布的各个区域之间移动，并且很快“忘记”它的历史状态。

• 无论是 $δ$ 太小（接受率高但步子小）还是 $δ$ 太大（接受率低总在原地踏步），相邻的样本都会非常相似，导致样本序列具有很高的自相关性。
• 高自相关意味着链的混合很差，我们需要生成更多的样本才能获得对目标分布的有效描述。
• 我们的目标: 是选择一个“恰到好处”的 $δ$，使得接受率不太高也不太低（经验上20%-50%之间较好），从而让链能够以最快的效率探索整个分布，即使得样本的自相关性最小化。这就是MCMC中的参数调优。

5.7. 独立Metropolis算法

这是Metropolis-Hastings算法的另一个非常重要的特例，该算法的核心特征是：它的提议分布 $q(x,y)$ 是一个只与目标位置 $y$ 有关的固定分布 $g(y)$。 $$q(x,y)=g(y)$$ 在独立Metropolis算法中，我们的下一步提议是从一个固定的“地图” $g(y)$ 上随机“空降”一个点，这个过程完全独立于我们当前所在的位置 $x$。

因此，一个好的独立Metropolis算法，要求我们选择的提议分布 $g(y)$ 本身就是对目标分布 $f(x)$ 的一个比较好的近似。由于提议分布的特殊形式，M-H算法的接受概率公式也相应地发生了变化。

根据独立提议的定义，我们有 $q(x,y)=g(y)$，同理，反向转移的概率就是 $q(y,x)=g(x)$。 $$\alpha(x,y)=\min\begin{pmatrix}f(y)g(x)\f(x)g(y)\end{pmatrix}$$ 如果我们定义重要性权重 (importance weight) 为 $w(z)=\frac{f(z)}{g(z)}$ ，那么接受概率公式可以被重写为： $$\alpha(x,y)=\min\left(\frac{w(y)}{w(x)},1\right)$$ 这个公式告诉我们，算法的决策逻辑是比较新提议点 $y$ 和当前点 $x$ 的重要性权重。如果新提议点的权重更高，那么我们就更倾向于接受这个移动。这就像一个MCMC版本的“人往高处走”，只不过这里的“高处”是由重要性权重来定义的。

该算法的完整流程。 从当前状态 $x_{n−1}​=x$ 开始：

1. 提议 (Proposal):
   • 从一个固定的提议分布 $g(y)$ 中抽取一个候选点 $y$。（注意：这一步与当前位置 $x$ 无关）
2. 决策 (Decision):
   • 生成一个从0到1的均匀分布随机数 $u∼Unif[0,1]$。
3. 更新 (Update):
   • 判断是否 $u\leq\min(\frac{f(y)g(x)}{f(x)g(y)},1)$。
   • 如果成立，则接受提议，令下一个状态 $x_n​=y$。
   • 如果不成立，则拒绝提议，令下一个状态保持不变 $x_n​=x$。

• 优点:
  • 如果我们能找到一个非常接近 $f(x)$ 的提议分布 $g(y)$，那么这个算法会极其高效。因为每次的提议都是“全局的”，而不是局部的微小移动，所以链可以很快地探索整个分布空间，并且样本之间的自相关性会非常低。
• 缺点:
  • 算法的性能非常依赖于我们选择的 $g(y)$ 的好坏。
  • 如果 $g(y)$ 对 $f(x)$ 的近似很差，特别是当 $g(y)$ 的尾部比 $f(x)$ 的尾部“更轻”（即下降得更快）时，算法性能会变得非常糟糕。链可能会偶然跳到一个 $f(x)$ 不小但 $g(x)$ 极小的点，导致该点的重要性权重 $w(x)$ 极大，从而使得链被“困”在这个点上很久都无法移动。

我们同样看一个例子：

我们的目标依然是生成服从标准正态分布的样本，其概率密度函数为 $ϕ(⋅)$。这次我们选择了一个更高级的提议分布：学生t分布 (Student's t distribution)，具体来说是自由度为3的t分布。

选择t分布的理由：

• 形状相似: t分布和正态分布一样，都是钟形曲线，关于0点对称。所以用t分布来近似正态分布，在形状上是合理的。
• 尾部更“重” (Heavier Tails): 这是最关键的一点！$t$分布的尾部比正态分布的尾部下降得更慢。这意味着在远离中心的地方，$t$分布的概率密度相对更大。这个“重尾”特性对于独立Metropolis算法至关重要，我们稍后会详细解释为什么。
• 提议: 从自由度为3的t分布 ($t_3$​) 中随机抽取一个候选点 $y$。
• 决策: 生成一个从0到1的均匀分布随机数 $u$。
• 更新: 使用独立Metropolis算法的接受概率公式进行判断，如果 $u\leq\min\left(\frac{\phi(y)t_3(x)}{\phi(x)t_3(y)},1\right)$ 则接受提议，令 $x_n​=y$；否则拒绝，令 $x_n​=x$。

一个非常重要的理论条件：独立Metropolis算法能够有效工作，当且仅当满足以下条件： $$C=\sup_x\frac{f(x)}{g(x)}<\infty$$ 这个公式的意思是，目标分布 $f(x)$ 与提议分布 $g(x)$ 的比率必须是有界的。也就是说，我们必须能找到一个常数$C$，使得对于所有的 $x$ ，$\frac{f(x)}{g(x)}$ 的值都不会超过它。

• 如果这个条件不满足，就意味着 $g(x)$ 的尾部比 $f(x)$ 的尾部“轻”太多（即 $g(x)$ 在尾部趋近于0的速度远快于 $f(x)$）。
• 这样可能会发生灾难性的情况：假设链偶然跳到了一个位于 $f(x)$ 尾部的点 $x_{stuck}​$，这个点的 $f(x_{stuck})$ 虽然小，但 $g(x_{stuck})$ 可能小到几乎为零。这会导致该点的重要性权重 $w(x_{stuck})=f(x_{stuck})/g(x_{stuck})$ 变得极其巨大。
• 之后，无论算法再怎么从 $g(y)$ 中提议新的候选点 $y$（这些新点很可能落在分布的中心区域），新点的重要性权重 $w(y)$ 都会远远小于那个巨大的 w$(x_{stuck}​)$。

结果就是，接受概率 $\alpha=\min(\frac{w(y)}{w(x_{stuck})},1)$ 会变得极小，导致算法几乎永远无法接受新的提议，链就会被“困”在 $x_{stuck}$​ 这个点上，彻底失去了探索整个分布的能力。

6. 吉布斯采样 (The Gibbs Sampler)

• 当我们面对一个多维联合分布 $f(x_1,\ldots,x_K)$ 时，直接从中抽样通常非常困难，特别是当这个分布形式复杂且非标准时。
• 我们的目标，就是为这个包含 $K$ 个随机变量的向量 $(X_1,\ldots,X_K)$ 生成样本。
• 变量间的相关性: 在多维空间中，各个变量之间往往存在着复杂的依赖和相关关系，使得我们无法将它们拆开独立处理。
• 高维空间的复杂性: 高维分布的几何形状可能非常奇怪，比如狭长的山脊、多个峰值（多模态）等。
• “维数灾难”: 如果我们尝试使用之前学的Metropolis-Hastings算法一次性对整个向量 $(X_1,\ldots,X_K)$ 进行提议和更新，那么这个“一步到位”的提议被接受的概率可能会随着维度 $K$ 的增加而急剧下降，导致算法效率极低。

吉布斯采样的精髓可以概括为“分而治之”。它没有像M-H算法那样对整个向量 $(X_1,\ldots,X_K)$ 进行“提议-决策”，而是将一个复杂的高维采样问题，分解成了一系列简单的一维采样问题。

具体做法是：轮流地、每次只对一个变量进行采样，并且在采样时，把所有其他变量都固定在当前值。

这个过程依赖于一个关键概念——全条件分布 (Full Conditional Distributions)。

6.1. 全条件分布

对于多维随机向量中的某一个变量 $X_k$​，它的全条件分布是指在给定所有其他变量 $X_j(j\neq k)$ 的取值时，这个变量 $X_k$​ 所服从的单变量条件分布。

其数学形式为： $$f(x_k|x_1,\ldots,x_{k-1},x_{k+1},\ldots,x_K),\quad k=1,\ldots,K$$ 果我们能够得到所有变量的全条件分布，那么那个复杂的多维联合分布采样问题，就可以被转化成一个依次从一系列简单的单变量条件分布中进行采样的序列过程。

在一些温和的条件下，这一整套的全条件分布，能够唯一地确定 那个原始的多维联合分布 $(X_1,\ldots,X_K)$，并因此也能唯一地确定所有单个变量的边缘分布。

6.2. 二维问题

设定是我们想要从一个二维的目标联合密度函数 $f_{xy}(x,y)$ 中采样，我们的随机向量是 $Z=(X,Y)$。

在这个二维案例中，我们有两个全条件分布：

• 给定 $x$ 时 $Y$ 的条件分布：$f_{y|x}(y|x)=\frac{f_{xy}(x,y)}{f_x(x)}$
• 给定 $y$ 时 $X$ 的条件分布：$f_{x|y}(x|y)=\frac{f_{xy}(x,y)}{f_y(y)}$

吉布斯采样能够运行的前提是，这两个一维的条件分布是我们已知的，并且可以很方便地从中进行抽样。

假设我们当前处在第 $n−1$ 步，状态为 $z^{(n-1)}=(x^{(n-1)},y^{(n-1)})$ 。要生成第 $n$ 步的状态 $z_n$​，我们需要执行以下两个步骤：

步骤 1: 更新 $x$

• 我们从 $X$ 的全条件分布中抽取一个新的样本 $x^{(n)}$。
• 在抽样时，我们将 $Y$ 的值固定在其前一步的取值 $y^{(n−1)}$ 上。 $$x^{(n)}\sim f_{x|y}(x|y^{(n-1)})$$ 步骤 2: 更新 $y$
• 我们从 $Y$ 的全条件分布中抽取一个新的样本 $y^{(n)}$。
• 关键点: 在抽样时，我们将 $X$ 的值固定在刚刚更新过的最新值 $x^{(n)}$ 上，而不是旧的值 $x^{(n−1)}$。 $$y^{(n)}\sim f_{y|x}(y|x^{(n)})$$ 这个步骤中，每一步都依赖于其他分量的最新值。这是吉布斯采样的一个核心特征。在更新 $y$ 时，我们没有使用旧的 $x^{(n−1)}$，而是立刻使用了在同一步骤中刚刚生成的 $x^{(n)}$。这种“即时更新”的策略是标准吉布斯采样的做法，它有助于加速算法的收敛。
• 一个“循环”包含了两步采样，但它构成了我们马尔可夫链的一次完整转移。
• 链的状态是二维向量 $z=(x,y)$。
• 从状态 $z_{(n−1)}$ 到 $z^{(n)}=(x^{(n)},y^{(n)})$ 的转移过程，其概率只依赖于 $z_{(n−1)}$，因此满足马尔可夫性质。
• 完整流程: 我们从一个任意的初始点 $z^{(0)}=(x^{(0)},y^{(0)})$ 开始，然后不断地重复上述的“两步循环”，生成一个状态序列 $z^{(0)},z^{(1)},z^{(2)},\ldots$。这个序列就是我们想要的马尔可夫链。在经过足够长的“burn-in”阶段后，这个序列中的样本就可以被看作是来自目标联合分布 $f_{xy}(x,y)$ 的样本。

6.3. K维情况

• 初始化: 算法首先需要一组初始值，为向量中的每一个变量都赋一个初始值：$x_1^{(0)},\ldots,x_K^{(0)}$。
• 一次完整的迭代 (从第0步到第1步): 这个过程包含 $K$ 个子步骤，依次对每一个变量进行更新。
  1. 更新第一个变量 $x_1$​:
     • 我们从 $x_1$​ 的全条件分布中抽样，生成新的值 $x_1^{(1)}$。
     • 在抽样时，所有其他变量都使用它们的旧值 (第0步的值)。
     • 公式为： $\mathrm{Sample~}x_1^{(1)}\mathrm{~from~}f(x_1|x_2^{(0)},x_3^{(0)},\ldots,x_K^{(0)})$。
  2. 更新第二个变量 $x_2$​:
     • 我们从 $x_2$​ 的全条件分布中抽样，生成新的值 $x_2^{(1)}$。
     • 关键点: 在为 $x_2$​ 抽样时，我们使用 $x_1$​ 刚刚更新过的最新值 $x_1^{(1)}$，而其他变量 $(x_3,\ldots,x_K)$ 则仍然使用它们的旧值。
     • 公式为： $\mathrm{Sample~}x_2^{(1)}\mathrm{~from~}f(x_2|x_1^{(1)},x_3^{(0)},\ldots,x_K^{(0)})$。
  3. ...... (以此类推)
  4. 更新最后一个变量 $x_K$​:
     • 当我们更新最后一个变量 $x_K​$ 时，所有在它之前的变量 $(x_1,\ldots,x_{K-1})$ 都已经在本轮迭代中被更新过了，所以我们全部使用它们的最新值。
     • 公式为：$\mathrm{Sample~}x_K^{(1)}\mathrm{~from~}f(x_K|x_1^{(1)},x_2^{(1)},\ldots,x_{K-1}^{(1)})$。
• 一次迭代: 完成上述从 $x_1​$ 到 $x_K$​ 的 $K$ 个步骤，就构成了一次完整的迭代。这次迭代的结果是，我们得到了一个新的样本向量 $(x_1^{(1)},\ldots,x_K^{(1)})$。这个新的向量就是马尔可夫链的下一个状态。
• 长期行为: 我们不断地重复这个包含K个步骤的迭代过程。经过大量的迭代之后，所生成的马尔可夫链（由一系列的样本向量构成），其联合分布会收敛到我们想要的目标联合分布 $f(X_1,\ldots,X_K)$。

6.4. MCMC实现与诊断

1. 可视化诊断工具：轨迹图 (Trace Plot)

• 样本路径图 (sample path)，也常被称为轨迹图 (trace plot) 或历史图 (history plot)。
• 它的画法很简单：横轴是迭代次数 $t$，纵轴是该次迭代生成的样本值 $x^{(t)}$。
  • “混合良好 (mixing well)”的轨迹图: 它看起来应该像一条“毛茸茸的毛毛虫”，在某个稳定的区间内剧烈地、随机地上下摆动，没有任何明显的趋势。这表明链已经收敛，并且正在有效地探索目标分布的高概率区域。
  • “混合很差 (mixing poorly)”的轨迹图: 它可能会呈现出缓慢的上升或下降趋势（说明还没收敛），或者长时间停留在一个值上，然后突然跳到另一个值，形成阶梯状。这说明链的移动非常困难，探索效率很低。

2. 预处理步骤：预烧/老化 (Burn-in)

• MCMC链是从一个任意的初始点开始的，它需要一段时间来“忘记”这个初始点，并逐渐游走到并进入平稳分布状态。这个早期的“游走”阶段的样本并不能代表我们的目标分布，因此必须丢弃。
• 我们观察轨迹图，确定链大约从第 D 次迭代开始进入了稳定的“毛毛虫”状态，然后我们就把前 D 个样本全部丢弃，只保留后续的样本用于分析和推断。

3.处理样本相关性：稀疏化 (Thinning)

• 问题: MCMC生成的样本是前后相关的，而不是独立同分布的。
• 策略: 稀疏化 (thinning) 的方法，即在burn-in之后，我们每隔 k-1 个样本再取一个（比如只保留第10、20、30...个样本），然后把中间的都丢弃掉。
  • 稀疏化的传统目的是为了降低样本的自相关性。
  • 然而，现在的统计学界普遍认为，为了进行参数估计（如计算后验均值），稀疏化是不必要的，甚至是有害的。因为它丢弃了大量有用的信息。虽然高相关性意味着有效样本量会减小，但保留所有样本通常比保留稀疏化后的子集能得到更精确的估计。
  • 现在，稀疏化更多是出于工程原因，比如当链非常长时，为了减少存储和计算的压力。

4.核心难题：多模态分布 (Multiple Modes)

• 问题: 这是单一MCMC链最难诊断的“陷阱”之一。如果目标分布有多个峰值（像驼峰一样的山脉），一条链很可能只在一个山峰（一个模式）附近探索，然后看起来完美地收敛了。但实际上，它完全错过了分布的其他重要部分。
• 解决方案：运行多条链 (Multiple Chains)
  • 从多个差异很大的初始值开始，独立地运行多条马尔可夫链。
  • 然后，我们比较链内 (within-chain) 的行为和链间 (between-chain) 的行为。
  • 如果收敛良好: 所有的链，无论从哪里开始，最终都应该收敛到同一个分布。它们的轨迹图会互相交织、重叠在一起，看起来就像是从同一个分布中抽取的。链间的方差应该和链内的方差差不多。
  • 如果收敛有问题: 如果有些链被困在了不同的模式中，那么它们的轨迹图会停留在完全不同的数值区间，永不相交。此时，链间的方差会远大于链内的方差。
  • 像Gelman-Rubin诊断（常被称为 R-hat ）这样的形式化诊断方法，其数学原理就是比较链内方差和链间方差，为我们提供一个量化的收敛判断标准。

![[Pasted image 20250927184429.png]]

![[Pasted image 20250927184449.png]]

• 上图 (理想的轨迹图): 这是一张混合良好 (well-mixed)、已经收敛的链的轨迹图。
  • 特征: 整个图形看起来像一条“毛茸茸的毛毛虫”，或者说一段没有明显趋势的白噪声。它在一个稳定的区间内（大约0.6到0.8）快速地上下波动。
  • 含义: 这表明链已经到达了它的平稳分布，并且正在有效地探索这个分布的核心区域。这是我们经过burn-in之后，希望看到的理想状态。
• 下图 (有问题的轨迹图): 这是一张非常典型的、混合很差 (poor-mixing) 并且带有漫长burn-in阶段的轨迹图。
  • 特征:
    1. 缓慢的趋势: 在最初的约5000次迭代中，链的值一直在缓慢地上升，有一个非常明显的大趋势。
    2. 高自相关: 整个路径看起来非常“平滑”，而不是像上图那样剧烈波动。这说明相邻的样本之间相关性非常高。
  • 含义: 这张图清晰地告诉我们，这条链从一个不太好的初始点（约0.3）开始，花了将近5000次迭代才慢慢“爬”到目标分布的中心区域。

6.5. 例子

6.5.1. 构建贝叶斯模型

1. 确定似然

我们观察到的数据 $y_1​,…,y_n$​ 被假设为独立同分布地来自于一个正态分布 $N(\mu,\tau^{-1})$ 。这里的 $μ$ 是未知的均值，$τ$ 是未知的精度 (precision)，即方差的倒数 $(\tau=1/\sigma^2)$ 。在贝叶斯模型中，使用精度而不是方差，常常会使数学推导更为简洁。

2. 选择先验

• 在贝叶斯框架下，我们需要为所有未知参数（这里是 $μ$ 和 $τ$）设定先验分布，以表达我们在看到数据之前的信念。
  • 对均值 μ，我们选择一个正态分布作为先验： $\mu\sim N(0,\omega^{-1})$。
  • 对精度 τ，我们选择一个伽马分布作为先验： $\tau\sim\mathrm{Gamma}(\alpha,\beta)$。
• 正态分布是正态似然（当方差已知时）的共轭先验，而伽马分布是正态似然（当均值已知时）的精度参数的共轭先验。

3. 目标与联合分布

• 我们的目标: 是从参数的联合后验分布 $p(\mu,\tau|\mathbf{y})$ 中进行抽样。这个后验分布综合了数据信息和我们的先验信念，包含了关于参数的所有知识。
• 联合分布: 模型中所有变量的完整联合分布$f_{y,\mu,\tau}(\mathbf{y},\mu,\tau)$，它等于似然 × $μ$的先验 × $τ$的先验。这是我们后续推导全条件分布的出发点。

6.5.2. 推导全条件分布

1. 推导 $τ$ 的全条件分布

• 方法: 我们把联合分布的表达式看作是变量 $τ$ 的函数，忽略所有不含 $τ$ 的项，只保留与 $τ$ 相关的部分。
• 结果: 整理后我们得到 $\pi_{\tau|\mu,y}(\tau|\mu,\mathbf{y})\propto\tau^{\alpha+\frac{n}{2}-1}\exp\left{-\left(\frac{\sum_{i=1}^n(y_i-\mu)^2}{2}+\beta\right)\tau\right}$
• 识别分布: 这个函数形式正是伽马分布的核。因此，我们可以直接写出 $τ$ 的全条件分布： $$\tau|\mu,\mathbf{y}\sim\mathrm{Gamma}\left(\alpha+\frac{n}{2},\frac{\sum_{i=1}^n(y_i-\mu)^2}{2}+\beta\right)$$ 2. 推导 $μ$ 的全条件分布
• 方法: 同样地，我们回到联合分布，这次只保留与 $μ$ 相关的项。
• 结果: 整理后我们得到 $\pi_{\mu|\tau,y}(\mu|\tau,\mathbf{y})\propto\exp\left{-\frac{\mu^2}{2}(n\tau+\omega)-\frac{2\mu(\sum y_i)\tau}{2}\right}$。
• 识别分布:
  • 看到指数上是变量 $μ$ 的二次函数，我们就应该立刻联想到正态分布。通过对指数部分进行配方法，我们可以将其整理成正态分布概率密度函数指数的标准形式，从而直接读出其均值和方差（或精度）。
  • 配方的结果，即 μ 的全条件分布是一个正态分布： $$\mu|\tau,\mathbf{y}\sim N\left(\bar{y}\frac{n\tau}{n\tau+\omega},\frac{1}{n\tau+\omega}\right)$$ 这个后验均值可以看作是数据信息（样本均值 $\overline{y}$​）和先验信息（先验均值0）的一个加权平均，其权重由数据精度 $nτ$ 和先验精度 $ω$ 决定。吉布斯采样器就在这两个全条件分布之间迭代。

6.5.3. 吉布斯采样循环

• 初始化 $\mu^{(0)},\tau^{(0)}$。
• 循环迭代 $t = 0, 1, 2, ...$
  • a. 从 $Gamma(…)$ 分布中抽样 $τ^{(t+1)}$，其中公式里的 $μ$ 使用当前值 $μ^{(t)}$。
  • b. 从 $N(…)$ 分布中抽样 $μ^{(t+1)}$，其中公式里的 $τ$ 使用刚刚更新的值 τ(t+1)。
• 经过burn-in后，收集到的 $(μ^{(t)},τ^{(t)})$ 对就是我们想要的来自联合后验分布的样本。

6.6. 先验

在前面的例子里，我们为模型参数 $μ$ 和 $τ$ 设定了先验分布的形式，但是这些先验分布自身的参数，即超参数 $ω,α,β$，应该如何设定？

核心原则：

• 除非我们对参数有非常强、非常可靠的先验知识，否则我们应该选择那些能使先验分布无信息性 的超参数。
• “无信息性”的先验也被称为模糊先验、平坦先验 或 扩散先验。

6.7. 吉布斯采样的有效性

证明我们通过吉布斯采样构建的马尔可夫链，其平稳分布 确实是我们想要的目标联合分布。

6.7.1. 设定问题和证明目标

1.定义转移核 (Transition Kernel):

• 首先，我们为二维吉布斯采样器定义转移核 $k(z^{(n-1)},z^{(n)})$。
• “转移核”是在连续状态空间中描述“转移概率”的工具，可以理解为从状态 $z^{(n−1)}$ 转移到状态 $z^{(n)}$ 的概率密度。
• 吉布斯采样的转移过程分两步：先从 $f_{x|y}(x|y^{(n-1)})$ 抽样得到 $x^{(n)}$，再从 $f_{y|x}(y|x^{(n)})$ 抽样得到 $y^{(n)}$。因此，完成从 $z^{(n−1)}$ 到 $z^{(n)}$ 这个完整转移的概率密度就是这两步概率密度的乘积： $$k(z^{(n-1)},z^{(n)})=f_{x|y}(x^{(n)}|y^{(n-1)})f_{y|x}(y^{(n)}|x^{(n)})$$ 2.设定证明目标:
• 要证明 $f_{xy}(x,y)$ 是该链的不变密度，我们需要证明它满足平稳分布的定义。在连续空间中，这个定义是一个积分方程：$\int k(z^{(n-1)},z^{(n)})f_{xy}(x^{(n-1)},y^{(n-1)})dx^{(n-1)}dy^{(n-1)}=f_{xy}(x^{(n)},y^{(n)})$
  • 左边: 假设当前系统已经处于平稳状态（即服从 $f_{xy}$​ 分布），我们从所有可能的旧状态 $z^{(n−1)}$ 出发，经过一步转移（乘以转移核 $k$），最终到达新状态 $z^{(n)}$ 的总概率密度。
  • 右边: 系统在平稳状态下，本身就处于新状态 $z^{(n)}$ 的概率密度。
  • 这个方程的含义是：对于一个平稳分布，流入任何一个状态的总概率，必须等于从这个状态流出的总概率，从而使该状态的概率保持不变。

6.7.2. 证明过程

$$\int\int f_{x|y}(x^{(n)}|y^{(n-1)})f_{y|x}(y^{(n)}|x^{(n)})f_{xy}(x^{(n-1)},y^{(n-1)})dx^{(n-1)}dy^{(n-1)}$$ 将 $f_{y|x}(y^{(n)}|x^{(n)})$ 移到积分外面，因为它与积分变量 $x^{(n-1)},y^{(n-1)}$ 无关 $$=f_{y|x}(y^{(n)}|x^{(n)})\int\int f_{x|y}(x^{(n)}|y^{(n-1)})f_{xy}(x^{(n-1)},y^{(n-1)})dx^{(n-1)}dy^{(n-1)}$$ 先对内层关于 $x^{(n-1)}$ 的积分进行处理。根据边缘分布的定义，$\int f_{xy}(x^{(n-1)},y^{(n-1)})dx^{(n-1)}=f_y(y^{(n-1)})$ $$=f_{y|x}(y^{(n)}|x^{(n)})\int f_{x|y}(x^{(n)}|y^{(n-1)})f_y(y^{(n-1)})dy^{(n-1)}$$ 使用条件概率的定义 $f_{x|y}\times f_y=f_{xy}$ ，将积分内的两项合并 $$=f_{y|x}(y^{(n)}|x^{(n)})\int f_{xy}(x^{(n)},y^{(n-1)})dy^{(n-1)}$$ 现在对 $y^{(n-1)}$ 进行积分。根据边缘分布的定义，这个积分的结果是 $f_x(x^{(n)})$ $$=f_{y|x}(y^{(n)}|x^{(n)})f_x(x^{(n)})$$ 再次使用条件概率的定义，我们知道 $f_{y|x}\times f_x=f_{xy}$ $$=f_{xy}(x^{(n)},y^{(n)})$$ 因此，二维吉布斯采样器拥有正确的不变密度（即 $f_{xy}$​）。只要我们构建的马尔可夫链是遍历的，那么链在第 $n$ 步的分布，将会收敛到这个目标联合分布 $f_{xy}​$。

6.8. 吉布斯采样是M-H算法的特例

在吉布斯采样的一个循环中，比如更新第 $i$ 个变量时，我们保持其他所有变量 $x_{−i}$​ 不变，只从全条件分布 $f_i(x_i|x_{-i})$ 中抽取一个新的值 $x_i^∗$​ 来更新第 $i$ 个分量。

们可以把这个过程看作是一次M-H提议。给定当前的状态向量 $X^{(t)}$，我们提议了一个候选向量 $X^∗$。这个候选向量 $X^∗$ 非常特殊：它与当前向量 $X^{(t)}$ 在所有其他分量上都完全相同，只在第 $i$ 个分量上可能不同。

为了让这个“伪装”成立，我们定义了一个巧妙的提议分布 $g_i​$。这个分布描述了我们如何抽取新的 $x_i^$ 。 $$g_i(x_i^|x_{-i}^{(t)})=\begin{cases}f_i(x_i^|x_{-i}^{(t)})&\text{如果 }x_{-i}^=x_{-i}^{(t)}\0&\text{其他情况}&\end{cases}$$

• 我们的提议方式是直接从第 $i$ 个变量的全条件分布 $f_i​$ 中进行抽样。
• 我们永远不会提议一个改变了其他任何分量的移动（因为在其他情况下，提议概率为0）。

代入M-H接受率公式，得到 $$\alpha=\frac{f(x^)/g_i(x_i^|x_{-i}^{(t)})}{f(x^{(t)})/g_i(x_i^{(t)}|x_{-i}^{(t)})}$$

任何联合分布 $f(x)$ 都可以被分解为条件分布和边缘分布的乘积，即 $f(x)=f_i(x_i|x_{-i})\times f(x_{-i})$ 。我们将这个洞察和 $g_i​$ 的定义代入接受率公式中： $$=\frac{f_i(x_i^|x_{-i}^{(t)})f(x_{-i}^{(t)})/f_i(x_i^|x_{-i}^{(t)})}{f_i(x_i^{(t)}|x_{-i}^{(t)})f(x_{-i}^{(t)})/f_i(x_i^{(t)}|x_{-i}^{(t)})}$$ 在分子和分母中，联合分布分解出的条件分布项，恰好可以和作为提议分布的条件分布项 $g_i​$ (也就是 $f_i$​) 完全抵消。 $$=\frac{f(x_{-i}^{(t)})}{f(x_{-i}^{(t)})}=1$$ 结果：接受概率 α 恒等于 1，所以每一个通过这种方式提议的候选点，都必然会被接受。这恰恰就是吉布斯采样的行为模式：从全条件分布中抽样，然后直接更新，没有任何“拒绝”的步骤。

这个证明告诉我们，吉布斯采样之所以不需要接受/拒绝步骤，是因为它所使用的提议分布（即全条件分布）是一种完美的提议分布，它使得M-H的接受率天然就达到了100%。因此，吉布斯采样是一种非常高效的MCMC方法，但它的应用前提是：我们必须能够直接从所有参数的全条件分布中进行抽样。

7. 数据增强

7.1. 核心思想

数据增强介绍的是一种策略，这个策略的核心是，当一个问题很难直接解决时，我们可以通过巧妙地引入一些“辅助”变量，让问题变得更容易处理。

数据增强的定义：

• 它是一种通过引入未观测数据 或潜变量 来构建迭代算法的方法。
• 它也被称为辅助变量法。

这个核心思想，在统计学界有两个著名的应用，分别对应着两种不同的目标：优化和采样。

A. 确定性算法 (Deterministic) - EM算法

• 在确定性算法中，这个思想被 Dempster, Laird, 和 Rubin (1977) 应用并推广，形成了 EM (Expectation-Maximization) 算法。
• 目标: EM算法的目标是优化，即寻找能让似然函数或后验概率密度达到最大值的参数（即求MLE或MAP估计）。
  • E-步 (Expectation): 这是“猜测”缺失部分的一步。在给定当前参数估计的情况下，计算这些潜变量的期望。
  • M-步 (Maximization): 这是“解决简单问题”的一步。我们假装上一步计算出的期望值就是潜变量的真实情况，将数据“补全”，然后在这个“完整数据”下轻松地计算出最大化似然的新参数估计。
  • 通过不断重复E-步和M-步，算法最终会收敛到参数的局部最优解。

B. 随机性算法 (Stochastic) - 数据增强算法

• 在随机性算法中，这个思想被 Tanner 和 Wong (1987) 推广，并直接命名为数据增强 (Data Augmentation) 算法。
• 目标: 这个算法的目标不是优化，而是采样 (Sampling)，即从一个复杂的目标后验分布中抽取样本。这正是MCMC的核心任务。
  • 数据增强算法可以被看作是EM算法的“随机”版本，它与吉布斯采样 (Gibbs Sampler) 密切相关。
  • I-步 (Imputation / 填充): 这是“随机猜测”缺失部分的一步。在给定当前参数样本的情况下，我们从潜变量的条件分布中随机抽取一个样本来“填充”数据。
  • P-步 (Posterior): 这是“随机解决简单问题”的一步。在给定刚刚填充好的“完整数据”下，我们从参数的条件后验分布中随机抽取一个新的参数样本（这一步通常很简单）。
  • 通过不断重复I-步和P-步，我们生成的参数样本序列就构成了一个马尔可夫链，其平稳分布就是我们想要的目标后验分布。实际上，这个过程就是一种将变量分为“参数”和“潜变量”两块的吉布斯采样。

7.2. 缺失数据机制

我们定义一系列记号，用来描述一个含有缺失值的数据集。

• $Y_{obs}​$: 我们能观测到的数据值。
• $Y_{mis}$​: 我们缺失的数据值。
• $Y=(Y_{obs​},Y_{mis}​)$: 这是一个概念上的“完整数据”，包含了观测部分和缺失部分。
• $θ$: 我们真正感兴趣的、描述数据生成过程的模型参数。

为了对“缺失”这个行为本身进行建模，我们定义一个缺失数据指示变量 $M$。

$M$ 是一个由0和1组成的向量或矩阵。如果对应的 $Y$ 值是观测到的，则 $M=1$；如果是缺失的，则 $M=0$。接下来我们可以把整个问题分解成两个部分：数据的生成过程和数据的缺失过程。

这个分解通过一个联合分布的因子分解来表示： $$f(Y,M|\theta,\psi)=f(Y|\theta)\times f(M|Y,\psi)$$

• $f(Y∣θ)$ - 数据生成模型:
  • 这是我们通常关心的科学模型部分。它描述了“完整数据” $Y$ 是如何从一个由参数 $θ$ 控制的分布中生成的。我们统计推断的目标，通常就是为了估计 $θ$。
• $f(M∣Y,ψ)$ - 缺失数据机制:
  • 这是描述“缺失”这个行为本身的模型。它告诉我们，一个数据点是缺失还是出现（由 $M$ 决定），其概率是否依赖于数据 $Y$ 的值。这个机制本身也可能由另一套我们不感兴趣的参数 $ψ$ 控制。

这个框架的重要性在于，它引导我们去思考一个关键问题：数据缺失的概率，到底依赖于什么？ 对这个问题的不同回答，引出了缺失数据三种著名的分类，也决定了我们后续应该如何处理数据：

1. 完全随机缺失 (Missing Completely at Random, MCAR):
   • 缺失的概率完全不依赖于数据 $Y$ 的任何值。即 $f(M∣Y,ψ)=f(M∣ψ)$。
   • 调查问卷的某一页因为打印机故障没印上，或者实验室的几个试管被意外打碎了。
2. 随机缺失 (Missing at Random, MAR):
   • 缺失的概率只依赖于观测到的数据 $Y_{obs}$​，而不依赖于缺失数据 $Y_{mis}​$ 本身。即 $f(M∣Y,ψ)=f(M∣Y_{obs}​,ψ)$。
   • 在一份健康调查中，男性比女性更不愿意回答关于心理健康的问题。这里，“是否回答问题”（$M$）依赖于“性别”（一个观测变量 $Y_{obs}​$），而不直接依赖于他真实的心理健康状况（可能是 $Y_{mis}​$）。
3. 非随机缺失 (Missing Not at Random, MNAR):
   • 缺失的概率依赖于缺失值 $Y_{mis}$​ 本身。
   • 例子：在收入调查中，收入非常高的人更倾向于不报告他们的收入。这里，“是否报告收入”（$M$）直接依赖于“收入本身”（$Y_{mis}​$）。这是最难处理的情况。

我们之所以要建立这个缺失数据机制的框架，就是为了判断我们的数据属于哪种情况。如果是MAR或MCAR，在某些条件下，我们可以使用一些比较简单的方法（比如数据增强）来处理，并且在分析时可以“忽略”掉缺失机制本身。但如果是MNAR，我们就必须对缺失机制 $f(M∣Y,ψ)$ 也进行建模，否则分析结果会有严重的偏差。

现在我们关心随机缺失的情况。首先，我们要知道如何从包含所有变量的联合分布，推导出我们实际拥有的数据（即观测数据 $Y_{obs}$​ 和缺失指示 $M$）的概率分布，这是通过对我们不知道的、缺失的数据 $Y_{mis}$​ 进行积分（或求和）来实现的，这个过程也叫作边缘化。 $$f(Y_{obs},M|\theta,\psi)=\int f(Y_{obs},Y_{mis}|\theta)\times f(M|Y_{obs},Y_{mis},\psi)dY_{mis}$$ 我们观测到某个结果的概率，等于所有可能导致这个结果的、我们未曾看到的“幕后故事”（即$Y_{mis}$的各种可能取值）的概率之和。在实际中，这个积分通常是非常复杂甚至无法计算的，这也是缺失数据问题的主要难点所在。

如果缺失数据的机制满足以下条件，我们就称数据是随机缺失的： $$f\left(M|Y_{obs},Y_{mis},\psi\right)=f(M|Y_{obs},\psi)$$ 当MAR假设成立时，在第一点的积分公式中，$f(M|Y_{obs},Y_{mis},\psi)$ 就可以被替换为 $f(M|Y_{obs},\psi)$ 。因为这一项与积分变量 $Y_{mis}​$ 无关，所以可以被提到积分号的外面。

整个观测数据的分布可以被完美地分解成两部分： $$f(Y_{obs},M|\theta,\psi)=f(M|Y_{obs},\psi)\times f(Y_{obs}|\theta)$$ 这个分解告诉我们，我们关心的参数 $θ$ 只出现在第二项 $f(Y_{obs}​∣θ)$ 中，而与缺失机制相关的参数 $ψ$ 只出现在第一项 $f(M∣Y_{obs}​,ψ)$ 中。为了对我们感兴趣的参数 $θ$ 进行推断（比如最大似然估计或贝叶斯推断），我们只需要关注 $f(Y_{obs}​∣θ)$ 这一部分就足够了，我们可以完全“忽略”掉那个描述数据为何会缺失的模型 $f(M∣Y_{obs}​,ψ)$。

因此，我们可以直接在观测到的数据上使用像最大似然法或数据增强这样的方法。这就是为什么MAR假设在现代统计学中如此核心和强大的原因。

正式的数学定义是，如果我们对主要参数 $θ$ 的贝叶斯推断，可以只通过观测数据的似然函数 $f(Y_{obs}​∣θ)$ 来进行，那么这个缺失数据的机制就是“可忽略的”。

要达到这种理想状态，需要满足两个条件：

1. 数据是随机缺失 (MAR): 它保证了观测数据的似然可以被完美地分解开。
2. 先验的独立性: 我们对主要参数 $θ$ 的先验信念，与我们对缺失机制参数 $ψ$ 的先验信念，必须是相互独立的。在数学上，这意味着联合先验可以分解为边缘先验的乘积：$π(θ,ψ)=π(θ)π(ψ)$。在绝大多数实际应用中，这个条件都是一个非常自然且合理的假设。

我们再来看完全随机缺失 (Missing Completely at Random, MCAR)。如果一个数据点是否缺失的概率，与所有的数据值（无论是观测到的 $Y_{obs}$​ 还是缺失的 $Y_{mis}​$）都完全无关，那么我们就称数据是完全随机缺失的。

数学公式: $$f(M|Y_{obs},Y_{mis},\psi)=f(M|\psi)$$ MCAR意味着数据缺失的发生是一个纯粹的随机事件，就像抛硬币一样，与数据本身的任何特质都无关。MCAR 是一个比 MAR 更强的假设。也就是说，如果数据是MCAR的，那么它必然是MAR的（因为如果缺失与任何数据都无关，那它自然也与缺失值本身无关）。反之则不成立。

由于 MCAR ⟹ MAR，所以如果数据是MCAR的（并且满足先验独立性），那么缺失机制一定是可忽略的。

• 当数据是MCAR时: 拥有完整数据的那些观测，可以看作是整个样本的一个完全随机的子集。在这种理想情况下，直接删除所有含有缺失值的行，然后只用剩下的完整数据进行分析，得到的结果在理论上是无偏的（尽管可能因为样本量减少而损失了一些效率）。
• 当数据是MAR但非MCAR时: 直接删除含有缺失值的行会导致偏差。比如，如果男性更多地不回答问题，我们删掉他们之后，样本里女性的比例就被人为地提高了，分析结果就会偏向女性的特征。此时，我们必须使用数据增强、吉布斯采样、EM算法或多重插补等更高级的方法来正确处理缺失值。

7.3. 数据增强 (DA) 算法

DA算法是EM算法的随机版本。

• EM算法的目标是优化，即寻找后验分布的峰值（MAP估计）或似然函数的最大值。它通过确定性的“期望-最大化”步骤，一步步“爬”到山顶。
• DA算法的目标是采样，即描绘出整个后验分布的全貌（而不仅仅是找到山顶）。它通过随机性的“填充-抽样”步骤，像一个MCMC算法一样在整座“山脉”上游走。

DA算法的基本思想：

• 策略性地引入潜变量 $Z$，使得两个关键的条件分布变得“可用” (available)。
• 这里的“可用”指的是，我们可以很方便地从这两个分布中进行抽样。
• 我们真正关心的是后验分布 $p(θ∣Y_{obs}​)$，但它通常因为含有复杂的积分而难以直接处理。
• 引入潜变量 $Z$ 的目的，就是把直接处理 $p(θ∣Y_{obs​})$ 这个“一步到位”的难题，分解成两个更容易处理的步骤。

DA算法的整个框架，都建立在这两个条件分布能够轻松采样的基础上。

1. 完整数据后验分布 $$f_{(\theta|Y_{obs},Z)}(\theta|Y_{obs},z)$$

• 这是参数 $θ$ 的后验分布，但它有一个重要的前提：假设我们已经知道了潜变量 $Z$ 的值。我们设计潜变量 $Z$ 的目的，就是为了让这个“补全了数据”之后的后验分布变得非常简单（比如变成一个标准正态分布或伽马分布，从而可以直接抽样）。

2. 条件预测分布 $$f_{(Z|Y_{obs},\theta)}(z|Y_{obs},\theta)$$

• 这是潜变量 $Z$ 的后验分布，其前提是假设我们已经知道了参数 $θ$ 的值。这个分布告诉我们，在当前的模型参数下，我们应该如何去“猜测”或“填充”那些缺失的信息 $Z$。
• DA算法就是一种吉布斯采样:
  • 这个框架实际上就是块吉布斯采样 (Blocked Gibbs Sampler) 的一个完美应用。
  • 我们将模型中所有未知量分成了两“块”：
    1. 参数块: $θ$
    2. 潜变量块: $Z$
  • DA算法的迭代过程，就是在这两个块的全条件分布之间轮流进行抽样。
    • 我们从 $f(Z∣Y_{obs}​,θ)$​ 中抽样 $Z$，这正是 $Z$ 的全条件分布 $p(Z∣θ,Y_{obs}​)$。
    • 我们从 $f(θ∣Y_{obs}​,Z)$​ 中抽样 $\theta$，这正是$\theta$的全条件分布 $p(θ∣Z,Y_{obs}​)$。
• 算法预览: 因此，DA算法的流程将会是一个包含两步的循环，不断迭代直至收敛：
  1. 填充步骤 (I-Step): 给定当前参数 $θ$ 的样本，从 $f(Z∣Y_{obs}​,θ)$​ 中抽取潜变量 $Z$ 的新样本。
  2. 后验步骤 (P-Step): 给定刚刚抽取的潜变量Z的新样本，从 $f(θ∣Y_{obs}​,Z)$​ 中抽取参数 $θ$ 的新样本。

这个迭代过程生成的 $(θ,Z)$ 序列，其平稳分布就是我们想求的联合后验分布 $p(θ,Z∣Y_{obs}​)$。我们只需保留其中的 $θ$ 样本，就完成了对目标后验 $p(θ∣Y_{obs}​)$ 的采样。

7.4. 原始DA算法

DA算法由两个交替进行的步骤组成：填充步骤 (imputation step / I-step) 和 后验步骤 (posterior step / P-step)。

1. 填充步骤 (I-step / Imputation Step)

• 这一步的任务是“猜测”或“填充”缺失的信息。但因为它是一个随机算法，所以它不是只猜一次，而是通过抽样生成多种可能的填充方案。
• 具体流程:
  1. 从我们当前对后验分布 $f_k(\theta|Y_{obs})$ 的估计中，抽取 $m$ 个参数样本 ${\theta^{(j)}}_{j=1}^m$​。
  2. 对于这 $m$ 个参数样本中的每一个 $θ^{(j)}$，我们都从潜变量的条件预测分布 $f_{(Z|Y_{obs},\theta)}(z|Y_{obs},\theta^{(j)})$ 中，抽取一个对应的潜变量样本 $z^{(j)}$。
• 这一步完成后，我们就得到了 $m$ 组“完整”的数据集，每一组都包含了观测数据 $Y_{obs}$​ 和一个我们“填充”进去的潜变量样本 $z^{(j)}$。
• 因为我们是通过多次（$m$次）抽样来填充缺失信息的，所以这个过程本身就构成了多重插补。多重插补的核心思想就是通过生成多个可能的完整数据集来恰当地反映由于数据缺失所带来的不确定性。

2. 后验步骤 (P-step / Posterior Step)

• 这一步的任务是利用我们刚刚“填充”好的 $m$ 组完整数据集，来更新我们对参数后验分布的估计。
• 更新公式:$$f_{k+1}(\theta|Y_{obs})=\frac{1}{m}\sum_{j=1}^mf_{(\theta|Y_{obs},Z)}(\theta|Y_{obs},z^{(j)})$$
  • 这个公式的含义非常直观。对于 m 组“完整”数据集中的每一组，我们都可以得到一个简单的“完整数据后验分布” $f_{(\theta|Y_{obs},Z)}(\theta|Y_{obs},z^{(j)})$。
  • 我们在第 $k+1$ 轮迭代中对后验分布的新估计，就是这 $m$ 个简单的后验分布的平均混合。我们通过对所有“可能的现实”（即 $m$ 种填充方案）进行平均，来得到一个更稳健、更准确的后验分布估计。

7.5. DA算法详情

核心目标: 假设我们的目标是计算形如 $∫f(x)p(x)dx$ 的积分，但这非常困难，因为我们无法直接从目标分布 $p(x)$ 中进行抽样。

• DA的策略:
  • 我们不直接处理这个棘手的边缘分布 $p(x)$，而是巧妙地引入一个辅助变量 $Y$，并构建一个联合分布 $p(x,y)$。
  • 这个联合分布 $p(x,y)$ 的构建需要满足两个条件：
    1. 当我们将辅助变量 $Y$ 通过积分“消除”掉之后，得到的边缘分布必须是我们最初的目标分布 $p(x)$，即 $∫p(x,y)dy=p(x)$。
    2. 这个联合分布的两个条件分布 $p(x∣y)$ 和 $p(y∣x)$ 必须是容易从中进行抽样的。
• 这就是数据增强思想的精髓。我们把一个难的边缘分布采样问题，转化成了一个简单的联合分布下的条件分布采样问题。

具体算法流程：

• 增强变量 (Augmented Variable):
  • 在这个框架下，我们引入的那个辅助变量 $Y$，就被称为增强变量。
  • 这个“增强”变量不一定具有真实的物理或现实意义。它很多时候纯粹是一个为了简化计算而引入的“数学工具”或“脚手架”。在我们之前讨论的DA算法版本中，“潜变量 $Z$”扮演的就是这个增强变量的角色。
• DA算法的两步循环:
  • 给定链的当前状态为 $x_n$​，算法通过以下两步来生成下一个状态 $x_{n+1}$​：
1. 第一步 (抽样$Y$): 固定 $X$ 的值为 $x_n​$，从条件分布 $p_{Y|X}(.|x_n)$ 中抽取一个增强变量 $Y$ 的样本，记为 $y_n$​。
2. 第二步 (抽样$X$): 固定 $Y$ 的值为刚刚抽出的 $y_n$​，从条件分布 $p_{X|Y}(.|y_n)$ 中抽取一个目标变量 $X$ 的样本，作为链的下一个状态 $X_{n+1}$​。
• 可以看到，整个系统有两个变量（或两“块”变量）：$X$ 和 $Y$。
• 算法的流程正是在这两个变量的全条件分布之间进行轮流抽样：
  • 第一步是从 $p(y∣x)$ 中抽样。
  • 第二步是从 $p(x∣y)$ 中抽样。

所以这个DA算法的核心和吉布斯采样是一样的。

7.6. 例子

我们的目标：从一个自由度为4的学生t分布 中进行抽样。其概率密度函数为： $$p(x)=\frac{3}{8}\left(1+\frac{x^2}{4}\right)^{-5/2}$$ 我们不直接处理 $p(x)$，而是引入一个辅助的增强变量 $y$，并精心构造一个二维的联合分布 $p(x,y)$： $$p(x,y)=\frac{4}{\sqrt{2\pi}}y^{3/2}\exp\left{-y\left(\frac{x^2}{2}+2\right)\right}$$ 这个看起来很复杂的 $p(x,y)$ 并非凭空捏造。它利用了 $t$ 分布的一个重要数学性质：$t$ 分布可以表示为正态分布的尺度混合。

简单来说，一个 $t$ 分布的随机变量，可以看作是一个正态分布的随机变量，除以一个伽马分布随机变量的平方根。我们引入的辅助变量 $y$ 在这里扮演的正是那个混合分布中精度 的角色。从这个构造出的联合分布 $p(x,y)$ 中，我们可以推导出两个非常简单的条件分布：

1. $X|Y=y\sim N(0,1/y)$ (给定y时，x服从一个正态分布)
2. $Y|X=x\sim\mathrm{Gamma}(5/2,x^2/2+2)$ (给定x时，y服从一个伽马分布)

• $p(x∣y)$**: 我们把联合分布 $p(x,y)$ 的表达式中所有不含 $x$ 的项都看作常数。与 $x$ 相关的项是 $\begin{aligned}\exp{-y\frac{x^2}{2}}=\exp{-\frac{x^2}{2(1/y)}}\end{aligned}$ 。这是均值为0、方差为 $1/y$ 的正态分布的核。
• $p(y∣x)$: 我们把 $p(x,y)$ 的表达式中所有不含 $y$ 的项都看作常数。与 $y$ 相关的项是 $y^{3/2}\exp{-y(\frac{x^2}{2}+2)}$ 。由于 $3/2=5/2−1$，这正是形状参数为 $5/2$、速率参数为 $(x^2/2+2)$ 的伽马分布的核。

通过数据增强，我们成功地将一个难的 $t$ 分布采样问题，转化成了一个简单的正态分布和伽马分布的轮流采样问题。

最终的吉布斯采样算法流程如下：

1. 随机初始化一个 $x_0$​ 的值。
2. 进入循环迭代 (对于 $n = 0, 1, 2, ...$):a. 第一步 (抽样$y$): 给定当前值 $x_n$​，从 $\mathrm{Gamma}(5/2,x_n^2/2+2)$ 分布中抽取一个新的 $y_n$​。 b. 第二步 (抽样$x$): 给定刚刚抽出的 $y_n$​，从 $N(0,1/y_n)$ 分布中抽取一个新的 $x_{n+1}​$。
3. 不断重复这个循环。在经过足够长的burn-in阶段后，我们收集到的样本序列 $x_1​,x_2​,x_3​,…$ 就是服从我们最初的目标自由度为4的t分布的样本。

7.7 分组计数数据

• 想象一个传感器，在连续的360个时间单位里，记录每个单位时间内通过的个体数量。
• 观测数据: 我们得到的观测结果是一个汇总表格：
  • 0 次通过: 出现了 139 次
  • 1 次通过: 出现了 128 次
  • 2 次通过: 出现了 55 次
  • 3 次通过: 出现了 25 次
  • 4次及以上通过 (4+): 出现了 13 次
• 问题就出在最后一项“4+”上。对于这13个观测，我们只知道其计数值大于等于4，但不知道确切的数值（可能是4, 5, 6, ...）。这种数据就是一种分组数据 或 截尾数据，可以看作是一种含有缺失信息的情况。

们假设在任何一个时间单位内，通过的个体数量服从一个泊松分布 $P(λ)$，其概率质量函数为： $$P(X=x)=\frac{e^{-\lambda}\lambda^x}{x!},\quad x=0,1,\ldots$$ 我们的目标就是根据观测到的数据，来估计这个未知的平均发生率 $λ$。所有观测事件发生概率的连乘积。我们有360个观测，分为5组。于计数值为0, 1, 2, 3的观测:

• 139个“0”的概率是 $(P(X=0))^{139}=(e^{-\lambda})^{139}$
• 128个“1”的概率是$(P(X=1))^{128}=(\frac{e^{-\lambda}\lambda^1}{1!})^{128}$
• 55个“2”的概率是$(P(X=2))^{55}=(\frac{e^{-\lambda}\lambda^2}{2!})^{55}$
• 25个“3”的概率是$(P(X=3))^{25}=(\frac{e^{-\lambda}\lambda^3}{3!})^{25}$

对于计数值为“4+”的观测，它的概率是 $P(X\geq4)=1-P(X<4)=1-(P(0)+P(1)+P(2)+P(3))$ ，也就是 $$1-\sum_{i=0}^3\frac{e^{-\lambda}\lambda^i}{i!}=1-e^{-\lambda}\sum_{i=0}^3\frac{\lambda^i}{i!}$$ 我们有13个这样的观测，所以这部分的贡献是 $(\ldots)^{13}$

把以上所有部分的概率乘起来，然后把与 $λ$ 相关的项进行合并，忽略掉与 $λ$ 无关的常数项（比如阶乘），就可以得到似然函数： $$L(\lambda|x_1,\ldots,x_5)\propto e^{-347\lambda}\lambda^{128+55\times2+25\times3}\left(1-e^{-\lambda}\sum_{i=0}^3\frac{\lambda^i}{i!}\right)^{13}$$ 这个观测数据的似然函数非常复杂，尤其是那个 $(1-\ldots)^{13}$ 项。直接对这个函数进行最大化来求 $λ$ 的最大似然估计（MLE）会很困难。在贝叶斯框架下，将它与一个先验相乘得到的后验分布会更加复杂，难以直接从中采样。

因此，我们的第一步是用“缺失”的观测值来“补全”数据。

• 这里的“缺失数据” Z，指的就是那13个计数值为“4+”的观测所对应的确切计数值。我们引入一个潜变量向量 $z=(z_1​,…,z_{13​})$ 来代表这13个未知的、但肯定大于等于4的整数。
• 这个操作就是数据增强的核心：通过引入潜变量，将一个不完整的数据问题，转化为一个更容易处理的完整数据问题。我们将数据分为两部分：
• $X$: 观测到的数据，即那347个计数值为0, 1, 2, 3的明确观测。
• $Z$: 缺失的数据，即那13个未知的确切计数值。

给定 $λ$ 时，$X$ 和 $Z$ 的分布:

• $X$ 的分布: 对于347个已知的观测值，它们的似然就是标准的泊松分布似然的连乘积： $$X|\lambda\sim\prod_{i=1}^{347}\frac{e^{-\lambda}\lambda^{x_i}}{x_i!}$$ Z的分布: 对于13个未知的计数值，我们知道它们都来自于同一个泊松分布，但有一个额外限制：它们的值必须大于等于4。因此，它们的分布是一个截断泊松分布。 $$Z|\lambda,x\sim\prod_{j=1}^{13}\frac{e^{-\lambda}\lambda^{z_j}}{z_j!}I(z_j\geq4)$$ 将X和Z的分布相乘，得到了“完整数据”的联合分布（或者叫“完整数据似然函数”）： $$\frac{e^{-360\lambda}\lambda^{\sum_ix_i+\sum_jz_j}}{\prod_ix_i!\prod_jz_j!}\prod_jI(z_j\geq4)$$ 注意观察这个“完整数据”的联合分布。如果我们暂时忽略掉与 λ 无关的常数项（分母），那么它作为 λ 的函数，其形式为 $e^{-360\lambda}\lambda^{\text{total count}}$。这正是一个标准的、包含360个观测值的泊松分布似然函数的形式。

我们现在为唯一的未知参数 $λ$ 设定一个先验分布，我们选择的先验是 $π(λ)∝1/λ$。这个先验被称为杰弗里斯先验 (Jeffreys prior)，是泊松分布速率参数 $λ$ 的一个标准无信息先验。选择它的目的是为了让先验知识对后验结果的影响降到最低。

从技术上讲，$π(λ)∝1/λ$ 等价于在 $log(λ)$ 上取一个平坦的先验，这对于一个大于0的速率参数来说是常见的做法。将这个先验与我们上一页得到的“完整数据似然函数”相乘，我们就得到了所有未知量（$λ$ 和 $Z$）的联合后验分布。 $$\frac{e^{-360\lambda}\lambda^{\sum_ix_i+\sum_jz_j}}{\prod_ix_i!\prod_jz_j!}\prod_jI(z_j\geq4)\times\frac{1}{\lambda}$$ 现在从复杂的联合后验分布出发，我们推导出每一个未知量（或每一块未知量）的全条件分布。

1. Z 的全条件分布 (缺失数据的分布)

• 结果: $Z|\lambda,x\sim\text{Truncated Poisson}(\lambda)$。
  • 这个结果非常直观。如果我们已经知道了事件发生的真实平均速率是 $λ$，那么对于那些我们只知道计数值 ≥4 的观测，它们的分布应该是什么？答案就是一个截断泊松分布 (Truncated Poisson)。也就是说，我们从一个标准的泊松分布 P(λ) 中抽样，但只保留那些值大于等于4的样本。
  • $Z$ 的条件分布实际上并不依赖于已知的观测值 $X$。这很自然，因为在泊松模型中，每个时间单位的计数值都是相互独立的，所以未知的计数值 $z_j​$ 只依赖于全局的发生率 $λ$，而与其他时间点的计数值 $x_i​$ 无关。

2. λ 的全条件分布 (模型参数的分布)

• 结果: $\lambda|x,z\sim\mathrm{Gamma}(\sum_ix_i+\sum_jz_j,360)$。
• 补充讲解 (这个Gamma分布是怎么来的？):
  • 我们知道，“完整数据”的似然函数是一个标准的泊松似然。而伽马分布 是泊松分布似然的共轭先验。
  • 即使我们使用的先验 $π(λ)∝1/λ$ 是一个不恰当的伽马分布（可以看作是 $Gamma(0,0)$），共轭性依然成立。
  • 共轭更新的规则非常简单，就是把先验的参数和数据的信息直接相加：
    • 后验形状参数 (shape) = 先验形状参数 + 总计数值 = $0+(\sum x_i+\sum z_j)$。
    • 后验速率参数 (rate) = 先验速率参数 + 总观测数 = $0+360$。
  • 这就直接得到了Gamma后验分布。已知数据的计数值总和：$\sum_ix_i=128\times1+55\times2+25\times3=313$

现在假设我们已经有了第 $t-1$ 次迭代的结果 $\lambda^{(t-1)}$

第 $t$ 次迭代包含以下两步：

1. I-步 (填充Z): 为13个缺失值生成新的样本 $Z_j^{(t)}$​。这些样本来自于以 $\lambda^{(t-1)}$ 为参数的截断泊松分布（$P(\lambda^{(t-1)})$，且值 ≥4）。
2. P-步 (更新$\lambda$): 利用刚刚生成的“完整数据”，从$\lambda$的后验分布（即一个伽马分布）中抽取新的样本 $λ(t)$： $$\mathrm{Simulate~}\lambda^{(t)}\sim\mathrm{Gamma}(313+\sum_{j=1}^{13}z_j^{(t)},360)$$ 在运行了成千上万次迭代（并去除了burn-in样本）后，我们得到了一系列的 $λ$ 样本 ${\lambda^{(1)},\lambda^{(2)},\ldots,\lambda^{(M)}}$。这些样本共同描绘了 λ 的后验分布。通常，我们会用这个后验分布的均值作为 λ 的点估计。PPT介绍了两种计算这个均值的方法。

方法一：经验平均值 (Empirical Average)

• 这是最直接、最符合直觉的方法。我们直接计算所有采集到的 $λ$ 样本的算术平均值。 $$\hat{\lambda}=\frac{1}{M}\sum_{i=1}^M\lambda^{(i)}$$ 这就是蒙特卡洛方法的核心思想：用样本的平均值来估计总体的期望值。

方法二：条件期望的平均值 (Rao-Blackwellization)

• 这是一种更高级、在统计上更优的方法。 $$\lambda=\frac{1}{M}\sum_{i=1}^ME[\lambda|x,z^{(i)}]=\frac{1}{360M}\sum_{i=1}^M(313+\sum_{j=1}^{13}z_j^{(i)})$$ 在每次迭代 $i$ 中，我们不直接使用那个随机抽取的、带有噪声的样本 $λ^{(i)}$，而是利用在这一步已知的其他信息（即 $z^{(i)}$）来计算 $λ$ 的条件期望 $E[λ∣x,z^{(i)}]$。然后，我们对这些条件期望值求平均。

在第 $i$ 次迭代的P-步中，我们知道 $λ$ 的全条件分布是一个 $\mathrm{Gamma}(313+\sum z_j^{(i)},360)$ 。而一个 $\mathrm{Gamma}(\alpha,\beta)$ 分布的期望值（均值）就是 $α/β$。因此，$E[\lambda|x,z^{(i)}]=\frac{313+\sum z_j^{(i)}}{360}$ 就是对这个期望值求平均。

这种使用条件期望来改进估计量的方法，被称为Rao-Blackwellization。Rao-Blackwell定理告诉我们，通过这种方式计算出的估计量 $\hat{λ}$，其方差会比直接使用经验平均值（方法一）得到的估计量更小。

直观的理解是，我们在每一步都用一个精确的、解析的期望值，来代替一个随机的、有波动的样本值。这就相当于在过程中“平滑”掉了一部分随机噪声，从而使得最终的估计结果更加稳定和精确。因此，只要参数的条件期望容易计算，方法二（Rao-Blackwellization）通常是比方法一更受青睐的选择。

8. Rao-Blackwellization

简单来说，Rao-Blackwellization是一种“免费”提升MCMC估计量精度的方法。我们的目标和标准的蒙特卡洛方法一样，都是为了估计某个期望值 $\mu=E_f[h(X)]$ 通常我们会生成一堆样本 $X_1,\ldots,X_n$ ，然后用经验平均值 $\hat{\mu}_{MC}=\frac{1}{n}\sum h(X_i)$ 来估计 $μ$。

Rao-Blackwellization的前提:

• 假设我们的样本 $X_i$​ 可以被分解为两个部分 $(X_{i1},X_{i2})$。
• 最关键的条件是，我们能够从数学上（解析地）计算出条件期望 $E[h(X_{i1})|X_{i2}]$。

Rao-Blackwellized估计量:

• 基于概率论中的全期望定律：$E[A]=E[E[A∣B]]$。
• 我们可以定义一个新的、经过改进的估计量，即Rao-Blackwellized估计量 $\hat{\mu}{RB}$： $$\hat{\mu}{RB}=\frac{1}{n}\sum_{i=1}^nE[h(X_{i1})|X_{i2}]$$
• 这个新估计量的核心思想是用“精确”代替“随机”。
• 在标准的 $\hat{\mu}_{MC}$ 中，我们平均的是一堆随机的函数值 $h(X_i​)$。
• 在 $\hat{\mu}{RB}$ 中，我们把计算分成了两部分，对于其中可以精确计算的部分（即条件期望 $E[h(X{i1})|X_{i2}]$），我们直接用其解析的、精确的数学期望来代替那个随机的、有噪声的样本值 $h(X_{i1})$。
• 我们相当于在最终求平均之前，预先“平滑”掉了一部分随机性。

$\hat{\mu}{RB}$ 和 $\hat{\mu}{MC}$ 具有相同的均值（都是无偏估计），但 $\hat{\mu}_{RB}$ - 具有更小的方差。方差越小，意味着估计量越稳定、越精确。

这个结论的证明基于概率论中的全方差定律： $$Var(A)=E[Var(A|B)]+Var(E[A|B])$$ 标准估计量的方差是 $$Var[\hat{\mu}{MC}]=\frac{1}{n}Var[h(X_i)]$$ 我们可以将 $Var[h(X_i)]$ 利用全方差定律分解： $$Var[h(X{i1})]=E[Var(h(X_{i1})|X_{i2})]+Var(E[h(X_{i1})|X_{i2})]$$ Rao-Blackwellized估计量的方差则是： $$Var[\hat{\mu}{RB}]=\frac{1}{n}Var(E[h(X{i1})|X_{i2}])$$ 所以有关系式： $$Var[\hat{\mu}{MC}]=\underbrace{\frac{1}{n}Var[E{h(X{1i})|X_{i2}}]}{Var[\hat{\mu}{RB}]}+\underbrace{\frac{1}{n}E[Var(h(X_{1i})|X_{i2})]}{\geq0}$$ 由于方差永远是非负的，所以 $E[Var(h(X{1i})|X_{i2})]$ 这一项必然大于等于0。因此，$Var[\hat{\mu}{MC}]\geq Var[\hat{\mu}{RB}]$ 。这就严格证明了Rao-Blackwellized估计量的方差更小。

数据增强（或吉布斯采样）框架下，Rao-Blackwellization是有天然的优势的。因为在数据增强的设定中，我们经常会遇到一种理想情况：参数 $θ$ 在给定“完整数据” $X_{com​}$（即观测数据+增强数据）下的后验分布 $θ∣X_{com}​$，是一个我们熟知的、具有明确的条件期望（均值）解析表达式的分布。

更深层的理论保证：充分统计量与UMVU估计

• 我们从统计理论层面，解释为什么这个通过条件期望得到的估计量如此之好。
• 核心结论: 如果“完整数据” $X_{com}​$ 包含了参数 $θ$ 的充分统计量 (sufficient statistic)，那么根据Rao-Blackwell定理，我们通过条件期望得到的这个估计量，就是该参数的UMVU估计量。
  • 充分统计量 (Sufficient Statistic):
    • 这是一个高度凝练了数据信息的统计量。如果一个统计量是“充分”的，意味着它包含了样本中关于未知参数 $θ$ 的全部信息。一旦我们知道了这个充分统计量的值，原始的、更庞杂的数据本身对于推断 $θ$ 就不再提供任何额外的信息。
    • 例子: 对于一个正态分布，样本均值和样本方差就是其参数（总体均值和方差）的充分统计量。
  • UMVU估计量 (UMVU estimator):
    • 这是一致最小方差无偏估计量 的缩写。它是统计估计中的“黄金标准”。
    • 无偏 (Unbiased): 指这个估计量的期望值等于参数的真实值。
    • 一致最小方差 (Uniformly Minimum-Variance): 指在所有无偏估计量的大家族中，UMVU估计量的方差是最小的。
    • 简单来说，UMVU估计量是在不引入系统性偏差的前提下，我们能得到的最精确、最稳定的估计量。

在数据增强的框架下使用Rao-Blackwellization，其意义远不止是“降低方差”这么简单。在很多设计良好的模型中（特别是使用共轭先验时），我们引入的增强数据 $X_{mis}​$ 与原始数据 $X_{obs}​$ 结合起来构成的“完整数据” $X_{com}​$，恰好就包含了我们想估计的参数 $θ$ 的充分统计量。

这种情况下，我们通过计算条件期望 $E[\theta|X_{com}]$ 得到的Rao-Blackwellized估计量，不仅仅是一个“更好”的估计量，它在理论上就是“最好”的无偏估计量。

注意：两个估计量都是无偏的，也就是当数据量很大的时候，两个数都能收敛到真正的平均值，但是后者的方差会显著的低。

9. 模拟退火 (Simulated Annealing)

9.1. 思想

我们的目标是找到一个函数的最大值，即找到使 $h(θ)$ 最大的那个 $θ$。一个常见的应用就是寻找最大似然估计 (MLE)。模拟退火的方法是构建并模拟一个马尔可夫链，这个链的平稳分布被特殊地设计为： $$\pi(\theta)\propto\exp(h(\theta)/T)$$ 这里的 $T$ 是一个我们引入的、可以控制的参数，被称为“温度 (temperature)”。

温度 $T$ 就像一个控制旋钮，通过改变它，我们可以调整算法在函数 $h$ 的“表面”上移动的快慢和自由度。

Metropolis算法的接受概率 $\alpha=\min(1,\frac{\pi(\theta_{new})}{\pi(\theta_{old})})$ 。代入我们模拟退火的平稳分布，这个接受概率就变成了 $\alpha=\min(1,\exp(\frac{h(\theta_{new})-h(\theta_{old})}{T}))$ 。

• 当T很高时（加热阶段）:
  • 分母 $T$ 很大，导致指数部分 $\frac{\Delta h}{T}$ 趋近于0，因此接受概率 $α$ 趋近于1。
  • 这意味着，几乎所有的移动都会被接受，无论是“上坡”（找到更好的解）还是“下坡”（接受一个更差的解）。
  • 这对应于金属被加热的状态，算法可以在整个解空间中自由地、随机地探索，从而有很大概率跳出局部最优解的陷阱。
• 当T很低时（冷却阶段）:
  • 分母 $T$ 很小。如果 $h(\theta_{new})<h(\theta_{old})$（一个“下坡”的移动），那么 $Δh$ 是一个负数，$\frac{\Delta h}{T}$ ​ 会是一个很大的负数，导致接受概率 $α$ 趋近于0。
  • 这意味着，算法此时几乎只接受“上坡”的移动，拒绝所有“下坡”的移动。它的行为变得像一个贪婪的爬山算法。
  • 随着 $T$ 逐渐趋近于0，平稳分布 $π(θ)$ 的概率会高度集中在函数 $h$ 的最大值所在的那个极小邻域内。

具体实现

这个实现本质上是一个经过特殊改造的Metropolis算法。

前提: 算法采用一个对称的转移链，即 $q(x,y)=q(y,x)$。这使得我们可以使用简化的Metropolis接受准则，而无需Hastings修正项。

• 初始化: 从一个随机的初始点 $θ_0$​ 开始。
• 提议: 在第 $n-1$ 步，我们通过在当前点 $θ_{n−1}​$ 的一个小邻域内进行随机游走，来生成一个候选点 $θ$。建议的方式是从一个均匀分布中抽样：$\theta\sim\mathrm{Unif}[\theta_{n-1}-\delta,\theta_{n-1}+\delta]$。
• 决策: 生成一个0到1之间的随机数 $u$。
• 接受/拒绝: 根据以下条件判断是否接受这个移动： $$u\leq\min\left(\frac{\exp(h(\theta)/T)}{\exp(h(\theta_{n-1})/T)},1\right)=\min(\exp(\Delta h/T),1)$$ 如果接受，则更新状态 $θ_n​=θ$；否则状态保持不变 $θ_n​=θ_{n−1}​$。

我们定义函数值的变化量为 $\Delta h=h(\theta)-h(\theta_{n-1})$

• 情况一：“上坡”（找到更好的解）:
  • 如果 $h(\theta)\geq h(\theta_{n-1})$，那么 $Δh≥0$。
  • 接受概率为 $\min(\exp(\text{非负数}/T),1)=1$。
  • 结论：任何能让目标函数值变好或保持不变的移动，都一定会被接受。
• 情况二：“下坡”（考虑一个更差的解）:
  • 如果 $h(\theta)<h(\theta_{n-1})$，那么 $Δh<0$。
  • 接受概率为 $\min(\exp(\text{负数}/T),1)=\exp(\Delta h/T)$，这是一个在0和1之间的数。
  • 结论：即使一个移动会使目标函数值变差，它也仍然有一定概率被接受。
• 核心作用: 正是这种有概率地接受“坏”解的机制，赋予了算法跳出局部最优解的能力。

一个完整的模拟退火算法框架，它包含一个外循环（降温）和一个内循环（Metropolis迭代）。

1. 初始化:
   • 设定一个足够高的初始温度 $T_{start}$​。
   • 随机生成一个初始解 $θ_0$​。
   • 设定一个极低的终止温度 $T_{stop}​$。
2. 外循环（降温过程）:
   • while 当前温度 $T>T_{stop}​$:
     • a. 内循环（在当前温度下达到热平衡）:
       • * for $i = 1$ to $L$ ($L$ 是每个温度下的迭代次数):
     • 执行Metropolis算法步骤（提议 -> 决策 -> 接受/拒绝）。
     • b. 降温: * 根据预设的退火策略，降低温度 $T$。例如，$T=c×T$，其中 $c$ 是小于1的常数（如0.99）。
3. 结束: 返回算法过程中找到的最优解 $θ$。