基础理论与核心思想

统计学的两大派系

“概率”是什么意思？两派答案不同——这会影响你怎么做推断、怎么解释区间和结论。想象一下，我们的目标是估计香港大学所有男生的平均身高。这个真实的平均身高，我们用 θ 表示。

频率派 (Frequentist)

频率派学者会这样想：

• 参数 θ 是唯一的、固定的：香港大学所有男生的平均身高是一个确定的数字（比如175.3cm），它就在那里，不会改变。它是“众神”视角下的真相，只是我们凡人不知道而已。
• 数据是随机的：我不可能测量每个人的身高，所以我随机抽取了100个男生作为样本。这个样本的平均身高（比如174.8cm）只是对真实值 θ 的一次估计。如果我重新抽100个人，会得到另一个略有不同的样本均值。所以，我的抽样过程是随机的。
• “95%置信区间”的真正含义：
  • 基于我这100人的样本，我算出了一个“95%置信区间”，比如说 [174.1cm, 175.5cm]。
  • 错误理解：认为“真实平均身高 θ 有95%的概率落在这个区间里”。
  • 正确理解 (频率派)：这个区间的构建过程（即抽样+计算）是随机的。这个“95%”描述的是我这个构建方法的长期成功率。如果全世界上有很多个像我一样的人，每个人都去港大独立抽100个男生，并用同样的方法计算一个置信区间，那么我们所有人计算出的区间中，大约有95%的区间会成功地“捕获”到那个唯一的、固定的真实值 θ。而我手里的这个区间，要么捕获到了，要么没捕获到，不存在“95%的概率”这一说，只是我更“相信”它捕获到了而已。

频率派认为：数据/统计量是随机的（因为会换样本），参数是固定的但未知。

贝叶斯派 (Bayesian)

贝叶斯学者会这样想：

• 参数 θ 是不确定的，所以它是随机变量：在我开始调查（抽样）之前，我不知道港大男生的平均身高 θ 到底是多少。它可能在170cm到180cm之间的任何地方，但我认为它在175cm左右的可能性更大。我的这种不确定性，可以用一个概率分布来描述，这就是先验分布 (Prior Distribution)。这完全是基于我主观的信念。
• 数据是证据：现在，我随机抽取了100个男生，得到了他们的身高数据。这些数据是新的“证据”。
• 更新信念：我使用贝叶斯定理，将我的先验信念 (Prior) 和新证据（数据似然 Likelihood） 结合起来，得到一个更新后的信念。这个更新后的信念也是一个概率分布，叫做后验分布 (Posterior Distribution)。
• “95%可信区间 (Credible Interval)”的含义：
  • 从后验分布中，我可以计算出一个“95%可信区间”，比如说 [174.2cm, 175.6cm]。
  • 这个区间的解释非常直观：“根据我已有的数据，真实平均身高 θ 有95%的可能性落在这个区间内。” 这正是大多数人直觉上对“置信区间”的错误理解。

贝叶斯派认为：把参数也当作随机变量，用先验表示不确定性。概率 = 主观信念的刻度。

贝叶斯推断的基石

贝叶斯定理：从先验到后验

核心公式：后验 ∝ 似然 × 先验

设 $A,BA,BA,B$ 是样本空间 $\Omega$ 里的两个事件，且 $P(B)>0$。 贝叶斯公式：$$P(A\mid B)=\frac{P(A\cap B)}{P(B)}=\frac{P(B\mid A)P(A)}{P(B)}$$

• $P(A)$：先验概率（在看到 $B$ 之前，你对 $A$ 的主观看法/频率信息）；
• $P(B\mid A)$：似然（如果 $A$ 真的发生，看到 $B$ 的可能性有多大）；
• $P(B)$：证据/边际概率（不管 $A$ 真不真，单看 $B$ 的整体概率）；
• $P(A\mid B)$：后验概率（在看到 $B$ 之后，你对 $A$ 的更新信念）。

如果 ${A_1,\ldots,A_m}$ 把样本空间“划分”为不相交且并起来为 $\Omega$ 的若干情形（叫分割/分区），那么对每个 $i$ 都有$$P(A_i\mid B)=\frac{P(B\mid A_i)P(A_i)}{\sum_{j=1}^mP(B\mid A_j)P(A_j)}$$ 分母就是全概率公式：$$P(B)=\sum_{j=1}^mP(B\mid A_j)P(A_j)$$ 把事件换成“参数 $\theta$”与“数据 $y$”：$$p(\theta\mid y)=\frac{p(y\mid\theta)p(\theta)}{p(y)},\quad p(y)=\int p(y\mid\theta)p(\theta)d\theta$$ 分母 $p(y)$ 叫边际似然/证据，难算时就需要后续课里讲的 Laplace 近似、重要性采样、MCMC 等数值方法来逼近。

四大要素：先验(Prior), 似然(Likelihood), 后验(Posterior), 证据(Evidence)

似然（Likelihood）和似然函数

假设我们观测到一组独立同分布（i.i.d.）的数据 $y={y_1,\dots,y_n}$ ，它们来自某个概率分布 $f(y\mid \theta)$，$\theta$ 是未知参数向量。这个分布由一个或多个我们不知道的参数 $θ$ 决定。

似然函数把“数据看作固定、参数看作变量”，定义是 $L(θ|y) = p(y|θ)$。由于数据是i.i.d.的，所以总的概率是每次观测概率的连乘积：$$L(\theta\mid y)\equiv p(y\mid\theta)=\prod_{i=1}^nf(y_i\mid\theta).$$ 它告诉我们：给定不同的 $\theta$ ，这些已经看到的数据有多“像”由它生成的。

似然函数问的是：“如果参数真的是 θ，那么我们观测到手头这组数据的可能性有多大？”它是一个关于 θ 的函数。这一点至关重要。在似然函数中，我们已经拿到了数据 y，所以 y 是固定的。我们要做的是尝试代入不同的 θ 值，看看哪个 θ 能让 $L(θ|y)$ 最大。

例子 (抛硬币)：假设我们抛了4次硬币，结果是 y = {正, 反, 正, 正}。

• 如果假设硬币是公平的 (θ = 0.5)，那么观测到这个结果的概率是 0.5 * 0.5 * 0.5 * 0.5 = 0.0625。
• 如果假设硬币有点偏，正面朝上的概率是 θ = 0.75，那么观测到这个结果的概率是 0.75 * (1-0.75) * 0.75 * 0.75 ≈ 0.105。
• 因为 0.105 > 0.0625，所以我们说，θ = 0.75 这个假设比 θ = 0.5 更有可能 (more likely) 解释我们的数据。

最大似然估计 (MLE)

最大似然估计(MLE)就是找到那个能让似然函数 L(θ|y) 取到最大值的 θ。同时，最大化 L(θ|y) 等价于最大化 log L(θ|y)（对数似然）。

在实际问题中，数据点通常是独立同分布的，总似然是每个数据点似然的乘积：$$L(\theta|X)=\prod_{i=1}^nP(x_i|\theta)$$ 这个连乘形式求导非常复杂。取对数后，就变成了加法： $$\ell(\theta|X)=\log(L(\theta|X))=\sum_{i=1}^n\log P(x_i|\theta)$$ 对连加形式求导就简单多了。同时，当数据集很大时，很多个小于1的概率值连乘，结果会是一个极小的数，计算机可能会因为它太接近于零而无法精确表示（数值下溢，underflow）。而对数似然是求和，数值上稳定得多。

在频率学派的角度来看，数据是随机的，参数是固定但未知的；所以我们实际是在挑一个在“长期重复抽样”里表现最好的点估计（$\hat\theta_{\text{MLE}}​$）。

先验分布与后验分布

贝叶斯派认为，我们已经观测到的数据 y 是固定的，它就是我们拥有的证据。而我们不确定的参数 θ被看作是随机的。那么根据贝叶斯公式，给定一个先验分布 p(θ)，参数 θ 的后验分布 p(θ|y) 由以下关系给出： $$p(\theta\mid y)\propto p(y\mid\theta)p(\theta).$$

• 这个公式就是用于更新信念的贝叶斯定理的精髓。
• p(θ|y) (后验)：我们最终的产出。在看到数据 y 之后，我们对 θ 的更新后的信念。
• p(y|θ) (似然)：这就是似然函数，贝叶斯派直接把它拿了过来，作为连接“参数假设”和“观测数据”的桥梁。它衡量了当前的 θ 对数据的解释程度。
• p(θ) (先验)：这是贝叶斯派独有的。它是在我们看到任何数据之前，对 θ 可能是什么样子的一个初始信念。比如，对于抛硬币，我们可能会假设 θ 在0.5附近的可能性比较大。
• ∝ (正比于)：因为 p(y|θ)p(θ) 计算出的结果不一定是一个合法的概率分布（积分不为1），所以我们先用正比符号。

为了让它成为一个合法的概率分布，我们需要除以一个归一化常数。这个常数就是 p(y|θ)p(θ) 在 θ 所有可能取值上的积分。 $$p(\theta|\mathbf{y})=\frac{p(\mathbf{y}|\theta)p(\theta)}{\int p(\mathbf{y}|\theta)p(\theta)d\theta}$$ 分母 ∫ p(y|θ)p(θ)dθ 就是我们前一页讲到的证据 p(y)。它代表了在我们的先验信念下，观测到当前这组数据的“边际可能性”有多大。

后验预测分布 (Posterior Predictive Distribution)

根据数据 y，来更新我们对参数 θ 的信念，得到后验分布 p(θ|y)。这属于参数推断（Inference）的范畴。而后验预测分布关心的是，在已经观测到数据 y 的条件下，对一个未来的新数据点 ỹ 的预测。

预测 ỹ 并非直接完成，而是通过参数 θ 作为桥梁。我们的预测综合了所有可能的参数 θ 值，并根据我们对这些 θ值的信任程度（即后验概率 p(θ|y)）对它们进行加权平均。最终的预测分布 f(ỹ|y) 是通过对 θ 进行积分（或“边缘化”）得到的： $$f(\tilde{y}|\mathbf{y})=\int f(\tilde{y}|\theta)p(\theta|\mathbf{y})d\theta$$ 继续使用之前抛硬币的例子。

假设我们想知道一枚硬币的真实偏倚（正面朝上的概率 θ）。

1. 推断阶段 (Inference)：我们先抛了10次 (y)，结果是“7正3反”。根据这些数据，我们更新了对 θ 的信念，得到了一个后验分布 p(θ|y)。这个分布可能告诉我们，θ 的值很可能在0.7附近，但它在0.6或0.8也有一定的可能性，而在0.1的可能性则非常小。
2. 预测阶段 (Prediction)：现在，我们想回答一个新问题：“下一次（第11次）抛硬币，结果是正面的概率有多大？” 这个问题就是关于未来数据 ỹ 的预测。

一个简单但不完备的方法是：

• 点估计法 (类似MLE)：我们从后验分布中找出一个最可能的 θ 值，比如 θ=0.7。然后直接用这个值做预测：“下次是正面的概率就是70%。”
• 这个方法的缺陷：它忽略了我们对 θ 的不确定性。我们的后验分布其实说的是 θ 可能是0.65，也可能是0.75，只是0.7的可能性最大而已。只用一个点，就丢失了大量信息。

贝叶斯派认为更严谨的做法是：

• 综合所有可能性：我们应该把所有可能的 θ 值都考虑进来。
  • 如果 θ 恰好是0.6，那么下次为正面的概率是60%。
  • 如果 θ 恰好是0.7，那么下次为正面的概率是70%。
  • 如果 θ 恰好是0.8，那么下次为正面的概率是80%。
  • ...等等
• 加权平均：我们该如何把这些预测融合起来呢？很简单，根据后验分布 p(θ|y) 给它们分配权重。后验分布告诉我们 θ=0.7 的可能性最大，θ=0.6 和 θ=0.8 的可能性次之。所以，我们最终的预测结果应该是这些单一预测的一个加权平均值。

所以再看这个公式： $$f(\tilde{y}|\mathbf{y})=\int f(\tilde{y}|\theta)p(\theta|\mathbf{y})d\theta$$

• f(ỹ|y)：我们想求的最终预测。在已知历史数据 y 的情况下，未来数据 ỹ 的概率分布。
• ∫ ... dθ：这个积分符号就代表了我们上面说的“对所有可能的θ进行加权求和/平均”的过程。
• f(ỹ|θ)：这是单一预测部分。如果我们知道参数 θ 的确切值，那么对未来数据 ỹ 的预测分布是什么。
• p(θ|y)：这是权重部分。也就是我们基于历史数据 y 得到的关于 θ 的后验分布，它告诉我们每个 θ 值的可信度有多高。

所以后验预测分布，就是用后验分布 p(θ|y) 作为权重，对由参数 θ 决定的所有单一预测 f(ỹ|θ) 进行的加权平均。这种方法不依赖任何一个单独的 θ 值，而是整合了所有 θ 的可能性。这使得贝叶斯的预测自然地包含了参数的不确定性。

贝叶斯推断实践

常用模型与共轭先验

共轭先验(Conjugate Prior)

假设我们有一组来自正态分布（钟形曲线）的数据，但我们不知道这个分布的均值 μ 是多少（比如，我们想知道用户的平均点击率，并假设它服从正态分布）。为了简化问题，我们假设分布的方差 σ² 是已知的。

我们经过以下三步贝叶斯处理：

• 选择先验 (Prior)：我们为未知的均值 μ 设定一个先验分布，PPT中选择的也是一个正态分布 N(μ₀, σ₀²)。这代表了我们在看到数据前对 μ 的初步猜测。

$$L(\mu|\mathbf{y})=\prod_{i=1}^n\frac{1}{\sqrt{2\pi}\sigma}\exp\left{-\frac{(y_i-\mu)^2}{2\sigma^2}\right}.$$

• 计算似然 (Likelihood)：根据我们收集到的数据 y，计算似然函数 L(μ|y)。
• 得到后验 (Posterior)：将“先验”和“似然”相乘，得到 μ 的后验分布 p(μ|y)。

$$\begin{aligned}p(\mu|\mathbf{y})&\propto\quad\exp\left{-\frac{\sum_{i=1}^n(y_i-\mu)^2}{2\sigma^2}\right}\times\exp\left{-\frac{(\mu-\mu_0)^2}{2\sigma_0^2}\right}\&\propto\quad\exp\left{-\frac{1}{2}{\left(\frac{\sigma_0^2+\sigma_n^2}{\sigma_0^2\sigma_n^2}\right)}{\left(\mu-\frac{\sigma_0^2\bar{y}_n+\mu_0\sigma_n^2}{\sigma_0^2+\sigma_n^2}\right)}^2\right},\end{aligned}$$

我们有如下结论：

• 当先验和似然都是正态分布时，最终的后验分布也仍然是一个正态分布。
• 这个新的后验分布的均值 $(\sigma_0^2\bar{y}_n+\mu_0\sigma_n^2)/(\sigma_0^2+\sigma_n^2)$ ，是先验均值 μ₀ 和数据样本均值 ȳₙ 的一个加权平均。方差是 $\sigma_0^2\sigma_n^2/(\sigma_0^2+\sigma_n^2)$
• 这个加权的权重，与各自的方差成反比。也就是说，谁的方差小（确定性高），谁的话语权就更大。

在这里，我们选择的先验是正态分布，而数据本身（似然函数的形式）也是源于正态分布。当这两者相乘后，得到的后验分布“恰好”还是一个正态分布。这种“先验分布”和“后验分布”属于同一个分布家族的优美特性，就叫做“共轭”。这种先验就叫共轭先验 (Conjugate Prior)。

这极大地简化了计算。我们不需要处理复杂的积分或者奇怪形状的后验分布，只需要根据公式，计算出新的正态分布的均值和方差就行了。

在这个例子中，后验均值的公式看起来很复杂，但它的本质就是一个加权平均： $$\text{后验均值}=w_\text{先验}\cdot\mu_0+w_\text{数据}\cdot\bar{y}_n$$

其中，权重 w 取决于我们对每一方信息的确定程度。在统计学中，方差的倒数（$1/σ^2$）常常被称为精度 (Precision)，它恰好可以用来衡量确定性。

• 先验的确定性：由先验方差 σ₀² 决定。σ₀² 越小，说明我们的先验知识越精确，我们对自己的初始猜测 μ₀就越自信。
• 数据的确定性：由数据样本均值的方差 σₙ² = σ²/n 决定。这个公式告诉我们两件事：
  1. 原始数据分布的方差 σ² 越小，数据的确定性越高。
  2. 数据量 n 越大，数据的确定性越高。这非常符合直觉：数据越多，真相越清晰。

Beta-Binomial 模型

二项分布 (Binomial Distribution)

• 二项分布 $Y ~ Bin(n, p)$ 描述的是：在 $n$ 次独立的伯努利试验中，事件“成功”发生的总次数。其中，每次试验成功的概率都是 $p$。
• 概率质量函数 (p.m.f.)：计算“在 $n$ 次试验中，恰好发生 $y$ 次成功”的概率的公式： $$P(Y=y)=\binom{n}{y}p^y(1-p)^{n-y}$$
• 二项分布的期望（平均值）和方差的计算公式：
  • 期望 $E(Y)=np$
  • 方差 $Var(Y)=np(1−p)$

贝塔分布 (Beta Distribution)

• 定义：贝塔分布 $X ~ Beta(α, β)$ 是一个定义在 (0, 1) 区间上的连续概率分布。它由两个正的形状参数 $α$ 和 $β$ 控制。$$f(x)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}x^{\alpha-1}(1-x)^{\beta-1},\quad0<x<1,$$
• 用途：它的值域是 (0, 1)，这使得它天生就适合用来为某个概率值 $p$ 进行建模。在贝叶斯统计中，它通常被用作描述概率参数的先验分布。
• 特殊情况：当 $α = 1$ 且 $β = 1$ 时，贝塔分布就退化为我们熟悉的 (0, 1) 上的均匀分布。这代表我们对 $p$ 的取值没有任何偏好，认为所有可能性都均等。
• 期望和方差：
  • 期望：$\mathsf{E}(X)=\frac{\alpha}{\alpha+\beta}$
  • 方差：$\mathsf{Var}(X)=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$

理解贝塔分布的关键是直观地理解两个形状参数 α 和 β 的作用。一个绝佳的理解方式是：

• α 可以看作是“伪成功次数” (pseudo-successes)
• β 可以看作是“伪失败次数” (pseudo-failures)

这两个参数共同描绘了我们在看到任何真实数据之前，心中对这个概率 p 的一个“假想”或“先验”的印象。

• α 和 β 的相对大小决定了分布的峰值位置（最可能的值）。
  • 如果 α > β，说明我们相信“成功”比“失败”多，分布的峰值会偏向1。
  • 如果 α < β，说明我们相信“失败”比“成功”多，分布的峰值会偏向0。
  • 期望公式 $\mathsf{E}(X)=\frac{\alpha}{\alpha+\beta}$ 完美地印证了这一点：期望值就是“成功次数”占“总次数”的比例。
• α 和 β 的总和 (α+β) 决定了分布的胖瘦（确定性程度）。
  • α+β 的值越大，分布就越“瘦高”，说明我们的先验信念非常笃定，我们对 p 的估计非常自信。
  • α+β 的值越小，分布就越“矮胖”，说明我们的信念很模糊，我们认为 p 的可能取值范围很广。

先验共轭

贝塔分布是二项分布的共轭先验。

• 先验 (Prior)：beta分布 $$\theta|y\thicksim\mathrm{Beta}(\alpha,\beta).$$
• 似然 (Likelihood)：二项分布 $$p(y|\theta)=\binom{n}{y}\theta^y(1-\theta)^{n-y}.$$
• 后验 (Posterior)：根据贝叶斯定理，后验 ∝ 似然 × 先验。一个Beta分布和一个Binomial分布的似然函数相乘后，得到的后验分布仍然是一个贝塔分布。 $$p(\theta|y)\propto p(y|\theta)p(\theta)\propto\theta^{\alpha+y-1}(1-\theta)^{\beta+n-y-1}.$$ $$\theta|y\thicksim\mathrm{Beta}(\alpha+y,\beta+n-y).$$

在实际应用中，先验可以根据实际来设置，因为有相信比率（期望）和总相信度（α + β），就可以量化到一个beta分布上。然后后验分布直接根据公式简单相加就可以。

Dirichlet-Multinomial 模型

多项分布似然

多项分布是之前讨论的二项分布 (Binomial Distribution) 的扩展和泛化，其描述了在 n 次独立的试验中，K 个互斥事件各自发生了多少次的概率。我们用用 pᵢ 表示事件 i 发生的概率，xᵢ表示事件 i 在 n 次试验中发生的次数。

那么，计算“在n次试验中，事件1发生x₁次，事件2发生x₂次，...，事件K发生xₖ次”这个特定组合的联合概率的公式如下： $$f(x_1,\ldots,x_K)=\frac{n!}{x_1!\cdots x_K!}p_1^{x_1}\cdots p_K^{x_K}$$ 如果我们只关心其中某一个事件 i 的发生次数 Xᵢ，那么它自己是服从二项分布 Bin(n, pᵢ) 的。

• 期望：$\mathsf{E}(X_i)=np_i$
• 方差：$\mathsf{Var}(X_i)=np_i(1-p_i)$
• 协方差：$\mathsf{Cov}(X_i,X_j)=-np_ip_j$

狄利克雷分布先验：贝塔分布的多维扩展

狄利克雷分布 (Dirichlet Distribution)可以看作是贝塔分布 (Beta Distribution) 的多维度扩展

• 定义：狄利克雷分布 X ~ Dir(α₁, ..., αₖ) 是一个向量的概率分布。这个向量 X = (x₁, ..., xₖ) 中的每个元素 xᵢ 都在0到1之间，并且所有元素加起来严格等于1。它由一组K个正的浓度参数 αᵢ 控制。
• 用途：它的定义完美匹配一个“K种互斥结果的概率集合”的性质（比如掷骰子时，6个面的概率），因此，它天生就是为多项分布的概率参数建模的理想先验分布。
• 与贝塔分布的关系：
  • 狄利克雷分布是贝塔分布的多变量版本。当 K=2 时，狄利克雷分布 Dir(α₁, α₂) 就等价于贝塔分布 Beta(α₁, α₂)。
  • 狄利克雷分布向量中的任何一个单独的元素 Xᵢ，其边际分布都服从一个贝塔分布 Beta(αᵢ, α - αᵢ)。

它定义在 概率向量 上，也就是 $(x_1​,x_2​,…,x_k​)$，其中每个 $x_i≥0$，并且 $$x_1+x_2+\cdots+x_k=1$$ 换句话说，Dirichlet 分布建模的对象是“在 $k$ 个结果之间分配概率”的分布。

记作： $$(x_1,x_2,\ldots,x_k)\sim\mathrm{Dir}(\alpha_1,\alpha_2,\ldots,\alpha_k)$$ 其中参数 $α_i>0$。

其概率密度函数 (PDF)为：$$f(x_1,\ldots,x_k;\alpha_1,\ldots,\alpha_k)=\frac{1}{B(\alpha)}\prod_{i=1}^kx_i^{\alpha_i-1},$$ 其中归一化常数是 多元 Beta 函数：$$B(\alpha)=\frac{\prod_{i=1}^k\Gamma(\alpha_i)}{\Gamma\left(\sum_{i=1}^k\alpha_i\right)},$$

• 每个 $α_i​$ 控制相应分量 $x_i$ 的“集中程度”。
• $α_i$ 越大，说明我们更“偏好”这个分量。
• $α=(α_1,…,α_k)$ 可以理解为 伪计数（pseudo-counts），像是先验观测次数。

记：$$\alpha_0=\sum_{i=1}^k\alpha_i.$$ 有期望：$$\mathbb{E}[X_i]=\frac{\alpha_i}{\alpha_0}.$$ 方差：$$\mathrm{Var}[X_i]=\frac{\alpha_i(\alpha_0-\alpha_i)}{\alpha_0^2(\alpha_0+1)}.$$

• 分子表示“某类 vs 其他类”的波动。
• 分母表示总的归一化程度。

协方差：$$\mathrm{Cov}[X_i,X_j]=-\frac{\alpha_i\alpha_j}{\alpha_0^2(\alpha_0+1)}.$$

后验更新：各类别的计数相加

• 先验 (Prior)：狄利克雷分布 $$\boldsymbol{\theta}|\mathbf{y}\sim\mathrm{Dir}(\alpha_1,\ldots,\alpha_K)$$
• 似然 (Likelihood)：多项分布。 $$L(\boldsymbol{\theta}|\mathbf{y})=\frac{n!}{y_1!\cdots y_K!}\theta_1^{y_1}\cdots\theta_K^{y_K}$$
• 后验 (Posterior)：由于共轭性，我们的后验信念仍然是一个狄利克雷分布 $$p(\boldsymbol{\theta}|\mathbf{y})\propto\theta_1^{\alpha_1+y_1-1}\cdots\theta_K^{\alpha_K+y_K-1}$$ $$\boldsymbol{\theta}|\mathbf{y}\sim\mathrm{Dir}(\alpha_1+y_1,\ldots,\alpha_K+y_K)$$

Normal-Normal 模型

• 先验 (Prior)：前文已述，我们给 $μ$ 一个正态先验：$\mu\sim\mathcal{N}(\mu_0,\tau^2)$
• 似然 (Likelihood)：假设观测数据$x_1,x_2,\ldots,x_n\sim\mathcal{N}(\mu,\sigma^2)$，似然函数是： $$p(x\mid\mu)=\prod_{i=1}^n\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x_i-\mu)^2}{2\sigma^2}\right)$$
• 后验 (Posterior)：由于共轭性，我们的后验信念仍然是一个正态分布。

其中参数更新公式为：

• 后验方差：$$\tau_n^2=\left(\frac{1}{\tau^2}+\frac{n}{\sigma^2}\right)^{-1}$$
• 后验均值：$$\mu_n=\tau_n^2\left(\frac{\mu_0}{\tau^2}+\frac{n\bar{x}}{\sigma^2}\right)$$ 后验均值 $μ_n$ 是 先验均值和样本均值的加权平均：$$\mu_n=w\mu_0+(1-w)\bar{x},\quad w=\frac{\sigma^2/n}{\tau^2+\sigma^2/n}$$ 也就是说，先验和数据会“折中”，谁更确定（方差更小）就更有话语权。

从后验分布到最终决策

后验分布 $f(θ|y)$ 是我们通过贝叶斯分析得到的“藏宝图”。它不是直接在地图上标一个“X”，而是告诉我们宝藏在每个位置的概率密度。有的地方概率高（山峰），有的地方概率低（平原）。这张图包含了我们关于宝藏位置（参数$θ$）的全部信息。

那么如果我们现在要做后验预测分布，我们要选择哪个点呢？

点估计

最大后验概率估计 (MAP)：寻找后验的“最高峰”

最大后验概率（MAP）估计，就是去后验分布这座“山脉”的最高峰（众数，mode）去寻找宝藏。这个点的θ值，是我们认为最有可能的那个值。

优点：

• 计算相对容易：寻找一个函数的最大值是一个优化问题，在很多高维的机器学习模型中，做优化比做积分（求均值）要容易得多。
• 与最大似然估计(MLE)的联系：这一点非常重要。我们这里可以用平坦先验 (Flat prior) π(θ) ∝ 1。
  • 平坦先验代表“我们对θ一无所知，认为所有θ的可能性都均等”。
  • 在这种情况下，贝叶斯公式变为：后验 ∝ 似然 × 1，即后验分布的形状完全由似然函数决定。
  • 因此，后验的最高峰（MAP）就等于似然的最高峰（MLE）。
  • 结论：频率派的MLE可以被看作是贝叶斯MAP在“无信息先验”下的一个特例。

后验均值：寻找后验的“重心”

给定数据 $x$，参数的后验分布是：$$p(\theta\mid x)=\frac{p(x\mid\theta)p(\theta)}{p(x)}$$ 那么 后验均值 就是这个分布的期望：$$\mathbb{E}[\theta\mid x]=\int\theta\operatorname{}p(\theta\mid x)\operatorname{}d\theta$$ 后验均值利用了整个后验分布的信息，而不仅仅是最高点。这一点很重要。在决策论中，如果我们把“估计误差”的损失函数定义为平方损失（即 Loss = (估计值 - 真实值)²），那么后验均值是能够使期望损失最小化的最优估计。这使得它在理论上具有非常好的性质。

对比

想象一下一个非对称的后验分布，比如个人收入分布，它通常是右偏的（少数极高收入者将平均值拉高）。

• 众数 (MAP)：会是人群中最常见的那个收入水平，可能并不高。
• 均值 (Posterior Mean)：会被那些极高的收入值“拉”向右边，数值上会比众数大。
• 如果后验分布是对称且单峰的（比如正态分布），那么均值、中位数和众数(MAP)都在同一个点，选择哪个都一样。
• 如果分布是偏斜的，或者你有特定的损失函数考量（比如你更不能容忍高估而不是低估），那么选择哪个就需要仔细考虑。
• 在很多复杂的机器学习模型中，由于后验均值（需要积分）计算困难，MAP因其计算上的便利性（优化问题）而被广泛采用。

计算统计工具

解析近似与蒙特卡洛模拟

解析近似方法

拉普拉斯近似 (Laplace Approximation)

在贝叶斯统计中，我们经常会遇到这样的后验分布公式：$$p(\theta|\mathbf{y})=\frac{p(\mathbf{y}|\theta)p(\theta)}{\int p(\mathbf{y}|\theta)p(\theta)d\theta}$$ 问题就出在分母的这个积分上。在大多数实际问题中，这个积分非常难算，甚至没有解析解。没有它，我们就无法得到一个标准化的后验分布。

核心思想：用正态分布拟合后验“山峰”

拉普拉斯近似的直觉如下：

1. 找到山峰：首先，我们找到我们想要近似的函数 f(x) 的最高点，也就是它的众数 (mode) x̂。
2. 测量山峰的形状：我们不仅关心山峰在哪里，还关心它是“尖峭”的还是“平缓”的。这个形状由函数在峰值点的曲率决定，数学上用二阶导数来衡量。
   • 二阶导数的绝对值越大，山峰越“尖”，说明函数在峰值附近下降得很快。
   • 二阶导数的绝对值越小，山峰越“缓”，说明函数比较平坦。
3. 构建完美的替代品：我们利用以上信息，构建一个正态分布（完美的钟形曲线）来替代原来的“小山峰”f(x)。
   • 正态分布的均值 μ 就设为原函数的峰值 x̂。
   • 正态分布的方差 σ² 就由原函数的曲率（二阶导数）来决定。山峰越尖，方差越小；山峰越缓，方差越大。具体来说，$\sigma^2=(-h^{\prime\prime}(\hat{x}))^{-1}$−1，其中 $h(x) = log(f(x))$。

现在，原来那个难算的对 f(x) 的积分，就变成了对一个简单的正态分布的积分，这个问题就迎刃而解了。

数学证明（低维）

假设我们的目标是近似计算下面这个积分：$$I=\int_{-\infty}^{\infty}f(z)dz$$ 在贝叶斯推断中，$f(z)$ 通常是一个未归一化的后验概率分布，比如 $f(z)=p(Data∣z)p(z)$。我们希望找到这个积分的值（即模型证据 $p(Data)）$，并用一个简单的分布来近似归一化后的$p(z)=\frac{f(z)}{I}$。

第一步：重写函数为指数形式

任何正函数 f(z) 都可以写成 $f(z)=exp(lnf(z))$。我们定义一个新的函数 $L(z)=lnf(z)$，也就是对数函数。那么，原积分变为：$$I=\int_{-\infty}^{\infty}\exp(L(z))dz$$ 第二步：找到函数的峰值

拉普拉斯近似的关键是在 $f(z)$ 的峰值（mode）周围进行近似。这个峰值对应于 $L(z)$ 的最大值点。我们把这个点记为 $z_0$​。要找到 $z_0$​，我们需要求解 $L(z)$ 的一阶导数为零的点：$$\left.\frac{dL(z)}{dz}\right|_{z=z_0}=0$$ 第三步：在峰值点进行二阶泰勒展开

我们在 $z_0$​ 附近对 $L(z)$ 进行泰勒展开，并取到二阶项：$$L(z)\approx L(z_0)+L^{\prime}(z_0)(z-z_0)+\frac{1}{2}L^{\prime\prime}(z_0)(z-z_0)^2$$

• 根据第二步的定义，$z_0$​ 是最大值点，所以一阶导数 $L′(z_0​)=0$。
• 因为 $z_0$​ 是一个最大值点，所以其二阶导数必须是负数，即 $L′′(z_0​)<0$。

为了方便，我们定义一个正数 $A=−L′′(z_0​)$。于是，泰勒展开简化为：$$L(z)\approx L(z_0)-\frac{1}{2}A(z-z_0)^2$$ 第四步：将展开式代入积分

现在，我们将这个近似的 $L(z)$ 代回到我们的积分 $I$ 中： $$I\approx\int_{-\infty}^{\infty}\exp\left(L(z_0)-\frac{1}{2}A(z-z_0)^2\right)dz$$ 我们可以把与积分变量 $z$ 无关的项 $exp(L(z_0​))$ 提出来：$$I\approx\exp(L(z_0))\int_{-\infty}^{\infty}\exp\left(-\frac{1}{2}A(z-z_0)^2\right)dz$$ 第五步：利用高斯积分求解

剩下的积分部分 $\int_{-\infty}^{\infty}\exp\left(-\frac{1}{2}A(z-z_0)^2\right)dz$ 是一个标准的高斯积分。

我们知道，一个高斯分布 $N(μ,σ2)$ 的概率密度函数 (PDF) 是： $$p(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$ 对比我们的积分项 $\exp\left(-\frac{1}{2}A(z-z_0)^2\right)$ 和高斯分布的指数部分，可以发现它们形式完全相同。通过匹配系数，我们得到：

• 均值 $μ=z_0​$
• $\frac{1}{\sigma^2}=A\mathrm{~}\Longrightarrow\mathrm{~}\sigma^2=\frac{1}{A}$

所以，我们的积分项可以看作是一个均值为 $z_0$​、方差为 $1/A$ 的高斯分布的核（没有归一化常数的部分）。因此： $$\int_{-\infty}^{\infty}\exp\left(-\frac{(z-z_0)^2}{2(1/A)}\right)dz=\sqrt{2\pi\sigma^2}=\sqrt{2\pi(1/A)}=\sqrt{\frac{2\pi}{A}}$$

第六步：得到最终近似结果

将第五步的结果代回到第四步的表达式中，我们得到积分 $I$ 的近似值： $$I\approx\exp(L(z_0))\sqrt{\frac{2\pi}{A}}=f(z_0)\sqrt{\frac{2\pi}{-L^{\prime\prime}(z_0)}}$$

这就是一维拉普拉斯近似对积分的求解结果。同时，我们也得到了对原概率分布 $p(z)=f(z)/I$ 的近似。因为$f(z)\approx\exp(L(z_0)-\frac{1}{2}A(z-z_0)^2)$ 这个形式正好是一个高斯分布，因此：$$p(z)\approx N(z|z_0,A^{-1})=N(z|z_0,[-L^{\prime\prime}(z_0)]^{-1})$$

应用：近似后验均值（MAP ≈ Mean）

在使用拉普拉斯近似的前提下，计算复杂的“后验均值”最终得到的结果，近似于我们更容易找到的“后验众数（MAP）

这其实不难理解，因为拉普拉斯近似的最终产出是一个高斯分布：$$p(\theta|D)\approx N(\theta|\theta_{MAP},H^{-1})$$ 因为高斯分布的均值和众数是同一个点，所以，当我们用这个高斯分布来近似原始的后验分布时，我们实际上已经做了一个核心假设：在这个近似的世界里，均值就等于众数。

这也揭示了拉普拉斯近似的一个核心特点或者说“代价”：它忽略了真实后验分布可能存在的偏斜（skewness）。

因为拉普拉斯近似的思想是：任何一个长得像“小山丘”的分布，我都可以用一个完美的正态分布去近似它。所以其一个结论就是后验均值 ≈ 后验众数 (MAP)。

这个等式成立的前提是拉普拉斯近似要足够好。这意味着，真实的后验分布本身就需要长得比较像一个对称的、单峰的正态分布。如果真实的后验分布是严重偏斜的，或者有多个峰值（多峰分布），那么它的“重心”（均值）和它的“最高峰”（众数/MAP）可能会相差很远，此时这个近似效果就会很差。

蒙特卡洛方法

核心思想：用抽样平均近似积分期望

Monte Carlo 方法的本质是 用随机样本来近似积分。比如我们想求期望：$$\mathbb{E}[f(\theta)]=\int f(\theta)p(\theta\mid D)d\theta$$ 如果能直接从 $p(\theta\mid D)$ 采样 $\theta^{(1)},\theta^{(2)},\ldots,\theta^{(N)}$ 那么有近似：$$\mathbb{E}[f(\theta)]\approx\frac{1}{N}\sum_{i=1}^Nf(\theta^{(i)})$$ 也就是说，用样本均值来代替积分。

理论基石：大数定律与中心极限定理

它的理论基石是大数定律 (Law of Large Numbers)。如果我们想计算一个函数 $h(x)$ 在$f(x)$ 分布下的期望值 $E{f​[h(x)]}$（本质上是一个积分），我们只需要从 $f(x)$ 中抽取大量随机样本 $x_1, ..., x_n$ ，然后计算这些样本 $h(x_i)$ 的算术平均值。当样本量 $n$ 足够大时，这个平均值就会收敛到真实的期望值。

蒙特卡洛求后验预测分布

如前文所述，通过对所有可能的参数 θ 进行积分（或加权平均）得到的：$$p(\tilde{y}|y)=\int p(\tilde{y}|\theta)p(\theta|y)d\theta$$ 假设已经通过一些方法（比如马尔可夫链蒙特卡洛 MCMC，如Metropolis-Hastings或Gibbs采样，或者变分推断 VI）从后验分布 p(θ∣y) 中获得了大量的样本。

我们称这些样本为：${\theta^{(1)},\theta^{(2)},\ldots,\theta^{(S)}}$ ，其中 $S$ 是样本的数量（比如 $S$ =10000）。

这里的每一个 $θ^{(s)}$ 都是一个从后验分布中抽出的、 plausible (貌似可信) 的参数值（或参数向量）。现在，计算后验预测分布的步骤如下：

第 1 步：从后验分布中抽取一个参数样本

从已经获得的后验样本集合中，随机抽取一个样本，我们称之为 $θ^{(s)}$。

第 2 步：使用这个参数样本生成一个新数据点

将 $θ^{(s)}$ 看作是模型的“真实”参数，然后从预测分布（也就是似然函数）中生成一个（或多个）新的数据点 $\tilde{y}^{(s)}$。 $$\tilde{y}^{(s)}\sim p(\tilde{y}|\theta^{(s)})$$ 第 3 步：重复以上过程

对所有的后验样本（从 $s$ = 1 到 $S$）重复第1步和第2步。这样，每一个后验参数样本 $θ^{(s)}$ 都会对应生成一个预测数据样本 $\tilde{y}^{(s)}$。

终结果

当你完成了这个过程，你就会得到一个包含 $S$ 个新数据点的集合： $${\tilde{y}^{(1)},\tilde{y}^{(2)},\ldots,\tilde{y}^{(S)}}$$ 这个集合，就是后验预测分布 $p(\tilde{y}|y)$ 的蒙特卡洛样本。这个样本集合非常有用，你可以用它来做各种分析：

• 可视化预测分布：
  • 画出 $\tilde{y}^{(s)}$ 的直方图 (Histogram) 或核密度估计 (Kernel Density Estimate, KDE) 图。这可以让你直观地看到新数据点最可能落在哪里，以及预测的不确定性有多大。
• 计算点估计：
  • 预测均值：计算样本的均值 $E[\tilde{y}|y]\approx\frac{1}{S}\sum_{s=1}^S\tilde{y}^{(s)}$。这可以作为对新数据的单点预测值。
  • 预测中位数：计算样本的中位数，这在分布偏斜时是比均值更稳健的估计。
• 构建预测区间 (Prediction Intervals)：
  • 这是最有用的功能之一，用来量化预测的不确定性。例如，要计算一个95%的预测区间，只需要找到样本的2.5%和97.5%分位数即可。这个区间告诉你，有95%的概率，下一个新的数据点会落在这个范围内。
• 后验预测检验 (Posterior Predictive Checks)：
  • 这是一种模型诊断的方法。你可以比较生成的预测样本 {$\tilde{y}^{(s)}$} 的分布与原始观测数据 $y$ 的分布。如果两者看起来很相似（例如，均值、方差、峰度等统计量接近），说明你的模型拟合得还不错。如果差异很大，则可能意味着你的模型假设有问题。

一个简单的例子

假设我们要预测下一次抛硬币的结果。

• 模型：伯努利分布，参数 $θ$ 是正面朝上的概率。
• 数据 ($y$)：我们观察了10次，其中有7次正面。
• 先验：我们假设一个$Beta(1, 1)$的均匀先验。
• 后验：根据贝叶斯定理，后验分布是 $p(θ∣y)=Beta(1+7,1+3)=Beta(8,4)$。

现在，我们用蒙特卡洛方法来得到下一次抛硬币结果 ($\tilde{y}$​) 的后验预测分布。

1. 从后验采样：从 $Beta(8,4)$ 分布中抽取大量的样本，比如10000个：${θ(1),θ(2),…,θ(10000)}$。
   • 例如，$θ(1)=0.68,θ(2)=0.75,θ(3)=0.61,…$
2. 生成新数据：对每一个 $θ(s)$，我们从伯努利分布 $Bernoulli(θ(s))$ 中生成一个新数据点 $\tilde{y}(s)$ (0代表反面，1代表正面)。
   • 对于 $θ(1)=0.68$，我们生成一个随机数，如果小于0.68，则 $\tilde{y}(1)=1$，否则为0。
   • 对于 $θ(2)=0.75$，我们生成一个随机数，如果小于0.75，则 $\tilde{y}​(2)=1$，否则为0。
   • ... 以此类推
3. 得到预测样本：最终我们得到一个由0和1组成的集合 ${\tilde{y}​(1),\tilde{y}​(2),…,\tilde{y}​(10000)}$。

分析结果： 我们可以计算这个集合中1的比例。这个比例就是下一次抛硬币结果为正面的后验预测概率。它约等于后验分布的均值 $E[θ∣y]=8/(8+4)=8/12≈0.67$。

蒙特卡洛积分的优点：高维度

传统的数值积分方法（如黎曼和、梯形法则）是在坐标轴上划分出均匀的网格，然后计算每个小格子的面积再相加。在一维上这很有效。但在二维上，网格点数是N²；三维是N³；一百维就是N¹⁰⁰。计算量会指数爆炸，这就是所谓“维度灾难”。

蒙特卡洛方法通过随机撒点，完全不受维度影响。它的收敛速度（误差减小的速度）与维度无关，始终是1/√n的级别。因此，对于高维积分，蒙特卡洛方法几乎是唯一可行的数值方法。

生成随机数

逆转换法 (Inversion Method)

CDF, $F(x)$：它回答的是“随机变量X取值小于等于x的概率是多少”，即  $$F(x) = P(X ≤ x)$$ 计算机通常只能生成最基本的[0, 1]区间上的均匀随机数，我们有从这种简单的随机数，生成服从正态分布、指数分布等更复杂分布的随机数的需求。

逆转换法（Inverse Transform / Inversion Method）是用均匀分布样本去“反推”任意目标分布样本的通用采样方法。核心是： 如果$U\sim\mathrm{Unif}(0,1)$，那么令 $$X=F^{-1}(U),$$ 就得到 $X$ 服从 CDF 为 $F$ 的分布（即目标分布）。

CDF 曲线在哪个区间比较“陡峭”，说明该区间的概率密度比较大。当我们从 $y$ 轴均匀抽样时，这些陡峭区间会“接住”更多的 $u$ 值，从而映射到更多的 $x$ 值，完美地复现了原始分布的概率密度。

理论上可以从任何分布中生成样本。但它有两个重要的前提：

1. 我们必须能够写出这个分布的累积分布函数F(x)的数学表达式。
2. 我们必须能够求出这个F(x)的反函数F⁻¹(u)的解析式。

转换法 (Transformation)

“转换法”不是指某一个单一、特定的方法，而是一大类通过对已知分布的随机变量进行函数变换，来生成目标分布的随机变量的技术的总称。

例子

伽马(Gamma) → 贝塔(Beta)

伽马函数，通常用大写希腊字母 $Γ$ 表示，它的标准定义是一个积分形式，由数学家欧拉提出： $$\Gamma(z)=\int_0^\infty t^{z-1}e^{-t}dt$$ 这个定义对于任何实部为正的复数 z 都成立。伽马函数最核心、最重要的一个性质是它的递推关系： $$\Gamma(z+1)=z\Gamma(z)$$ 这个性质看起来和阶乘的性质 $n!=n×(n−1)!$ 非常相似，正是通过这个性质，伽马函数和阶乘联系了起来。

而伽马分布是一个描述等待特定数量事件发生所需要的时间的连续概率分布。其密度函数为： $$f(x)=\frac{1}{\Gamma(\alpha)\theta^\alpha}x^{\alpha-1}e^{-x/\theta},\quad x>0$$

若把“每次发生一个小事件所需时间”看成指数分布，那么**$α$ 次事件的总等待时间**就是 $Gamma(α,θ)$。所以 $Gamma$ 可以理解为“很多个独立指数的和”。

参数：

• $α>0$：形状（shape）
• $θ>0$：尺度（scale）；或用 $λ=1/θ$ 叫做“率”（rate）

那么有结论：

若 $X\sim\mathrm{Gamma}(\alpha,\theta)\cdot Y\sim\mathrm{Gamma}(\beta,\theta)$ 相互独立，并且尺度相同（同一 $θ$，或用 rate 形式就是同一 $λ$），那么 $$U=\frac{X}{X+Y}\sim\mathrm{Beta}(\alpha,\beta)$$ 同时，总和 $V=X+Y\sim\mathrm{Gamma}(\alpha+\beta,\theta)$ ，而且 $U$ 和 $V$ 独立。

(伽马 + 泊松) → 负二项分布

泊松分布专门用来模拟在固定时间或空间内，某个事件发生的次数 $$\Pr(N=n\mid\lambda)=\frac{e^{-\lambda}\lambda^n}{n!},\quad n=0,1,2,\ldots$$ 其参数为速率 $\lambda$ ，性质为：$\mathbb{E}[N]=\mathrm{Var}(N)=\lambda\text{(均值=方差)}$

但是在真实世界中，数据的方差往往比均值大得多。比如，客服中心可能大部分小时接到8-12个电话，但如果公司产品出了问题，某个小时的电话量可能飙升到100个。这一个极端值会极大地拉高整体的方差，导致方差远远大于均值。这种现象就叫做“过度离散” (Overdispersion)。

负二项分布也是一个“计数”分布，和泊松分布是“亲戚”，都可以用来模拟事件发生的次数。

但是其有两个参数，而不是一个。这使得它可以将均值和方差分开，特别擅长处理方差大于均值（过度离散）的数据。可以把它看作是泊松分布的一个更强大、更通用的升级版。

$$\Pr(Y=n)=\binom{n+r-1}{n}(1-p)^rp^n,\quad n=0,1,2,\ldots$$ 其中 $0<p<1,\mathrm{~}r>0$

负二项分布最初的定义来自于这样一个场景：不断重复做一件事（比如抛硬币），这件事只有“成功”和“失败”两种结果。负二项分布描述的是：为了达到 $r$ 次成功，你需要经历多少次失败。 例如，“我需要抛多少次反面（失败），才能最终看到第5次正面（成功）？”

结论是先从伽马分布抽一个泊松强度，再用这个强度抽泊松计数： $$\lambda\sim\mathrm{Gamma}(r,\beta),\quad Y\mid\lambda\sim\mathrm{Poisson}(\lambda).$$ 则边缘分布 $Y$ 服从负二项： $$Y\sim\mathrm{NegBin}\left(r,p=\frac{1}{\beta+1}\right).$$

标准正态 → 多元正态

多元正态分布（Multivariate Normal Distribution, MVN）是正态分布这个钟形曲线在多个维度上的推广。它不再是描述一个变量，而是同时描述多个相互关联的随机变量的联合分布。

• 一维正态分布: 描述一个变量。例如，只看身高的分布。它由两个参数定义：
  • 均值 $μ$: 钟形曲线的中心位置。
  • 方差 $σ^2$: 钟形曲线的“胖瘦”程度（离散程度）。
• 多元正态分布: 描述多个变量。例如，同时看身高、体重、年龄的分布。它也由两个参数定义，但它们从数值升级为了向量和矩阵：
  • 均值向量 $μ$​: 这是分布的中心点。如果描述d个变量，它就是一个d维向量，每个元素是对应变量的均值。例如 [平均身高, 平均体重, 平均年龄]。
  • 协方差矩阵 $Σ$ (Sigma): 这是多元正态分布的灵魂，也是最关键的概念。它同时描述了每个变量自身的离散程度以及变量之间的相互关系。

协方差矩阵是一个 $d×d$ 的对称矩阵（$d$ 是变量的个数）。我们用一个二维的例子（比如身高和体重）来理解它：

假设我们有两个变量 $X_1$​ (身高) 和 $X_2$​ (体重)，协方差矩阵 $Σ$ 就是一个 2×2 的矩阵： $$\Sigma=\begin{pmatrix}\operatorname{Var}(X_1)&\operatorname{Cov}(X_1,X_2)\\operatorname{Cov}(X_2,X_1)&\operatorname{Var}(X_2)\end{pmatrix}$$

• 对角线元素 (Diagonal):
  • Var($X_1$​) 和 Var($X_2$​) 分别是身高和体重的方差。
  • 它们告诉你每个变量自己作用的时候有多分散。方差越大，这个变量的分布就越“胖”。
• 非对角线元素 (Off-diagonal):
  • $Cov(X_1​,X_2​)$ 是身高和体重的协方差。
  • 它描述了两个变量之间的线性关系。
    • 协方差 > 0: 表示正相关。一个变量增大时，另一个也倾向于增大。
    • 协方差 < 0: 表示负相关。一个变量增大时，另一个倾向于减小。
    • 协方差 = 0: 表示两者线性无关。

MVN的概率密度函数 (PDF)公式看起来有点复杂，但结构和一维的非常像： $$f(x;\mu,\Sigma)=\frac{1}{\sqrt{(2\pi)^d|\Sigma|}}\exp\left(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)\right)$$ 我们的目标是$\mathbf{X}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ ，而我们现在有的 $Z$ 是一个向量，它的每个元素都是独立的 $N(0, 1)$ 随机数。在几何上，$Z$ 对应的概率云是一个完美的、中心在原点的、没有被拉伸或旋转的“圆形”云。

在一维中，公式为：$x = μ + σz$ ，在高维空间，这个公式变成了：$$X=\mu+L\cdot Z$$ 其中 L 是一个起到了“拉伸和旋转”作用的矩阵，类似于一维里的 σ。我们知道，变换后的协方差矩阵 Σ 是由变换矩阵 L 决定的。数学上可以证明，新向量 LZ 的协方差矩阵是 L Lᵀ (L 乘以 L的转置)。 所以，我们想让最终的协方差是 Σ，就必须找到一个矩阵 L，使得： $$\Sigma=LL^T$$ 求这个矩阵常用Cholesky 分解，粗略地看，我们可以认为 L 是协方差矩阵 Σ 的“平方根”，就像 σ 是方差 σ² 的平方根一样。

经过这两步操作，向量 X 就成了一个完美的、服从我们目标分布 N(μ, Σ) 的随机向量。

序贯模拟 (Sequential Simulation)

假如我们的目标是从一个高维联合分布 $$P(X_1,X_2,\ldots,X_n)$$ 中抽取样本。这里的核心挑战是，变量之间通常是相互依赖、不独立的，这导致在很多情况下，直接从联合分布采样很困难，因为它可能没有简单的闭式形式，或者维度太高。

但如果这个联合分布可以分解为条件分布的链式展开： $$P(X_1,X_2,\ldots,X_n)=P(X_1)\cdot P(X_2\mid X_1)\cdot P(X_3\mid X_1,X_2)\cdots P(X_n\mid X_1,\ldots,X_{n-1}),$$ 那么我们就可以 依次生成每个变量，这就是序贯模拟的核心。

• 先从边缘分布 $f_x(x)$ 中抽取一个 $x_i$ 。
• 以第一步抽到的 $x_i$ 为条件，从条件分布  $f_{y|x}(y|x=x_i)$  中抽取一个 $y_i$ 。
• 再以 $x_i$ 和 $y_i$ 为条件，从条件分布  $f_{z|x,y}(z|x=x_i, y=y_i)$  中抽取一个 $z_i$。
• 最终得到的 $(x_i, y_i, z_i)$  就是联合分布  $f(x, y, z)$  的一个有效样本。

假设我们想生成一组模拟数据，描述某城市一天的状况，包含三个相关的变量：

• x: 早晨是否下雨 (是/否)
• y: 下午交通是否拥堵 (是/否)
• z: 晚上能否准时赴约 (是/否)

这三个变量显然不是独立的。下雨会影响交通，交通会影响赴约。我们想从它们的联合分布 P(x, y, z) 中抽样。

序贯模拟就是按照事件发生的自然顺序，一步一步来：

第一步：模拟早晨的天气 x

• 我们需要知道边缘概率 P(x)，即“在没有任何其他信息的情况下，早晨下雨的概率是多少？”
• 假设根据历史数据，P(x="下雨") = 0.3。
• 我们生成一个[0,1]的随机数，如果小于0.3，我们就设定 xᵢ = "下雨"。假设我们这次抽到了“下雨”。

第二步：模拟下午的交通 y

• 交通状况取决于早晨的天气。所以我们需要条件概率 P(y|x)。
• 比如，P(y="拥堵" | x="下雨") = 0.8，而 P(y="拥堵" | x="晴天") = 0.2。
• 因为我们第一步抽到了xᵢ="下雨"，所以我们现在基于80%的概率来抽样。生成一个随机数，如果小于0.8，我们就设定 yᵢ = "拥堵"。假设我们这次抽到了**“拥堵”**。

第三步：模拟晚上的约会 z

• 约会能否准时，可能同时取决于天气和交通。所以我们需要条件概率 P(z|x, y)。
• 比如，P(z="迟到" | x="下雨", y="拥堵") = 0.9。
• 因为我们前两步的结果是“下雨”和“拥堵”，我们现在就基于90%的概率来抽样，很可能得到 zᵢ = "迟到"。

最终样本：经过这三步，我们就得到了一个完整的、符合逻辑链条的样本：("下雨", "拥堵", "迟到")。重复这个过程一万次，我们就能得到一个非常逼真的模拟数据集。

序贯模拟的美妙之处在于它的直观性，但它的实用性有两大前提：

1. 你需要知道所有这些边缘和条件分布的解析形式。在上面的天气例子和所有贝叶斯网络中，模型本身就是用这些条件概率定义的，所以这个前提成立。但在很多其他模型（比如物理系统、复杂的经济模型）中，我们可能只知道联合分布f(x,y,z)的数学形式（甚至只知道它正比于什么），但要从中解析地推导出f(x)和f(y|x)等可能会非常困难，甚至不可能。
2. 你需要能够从这些条件分布中进行抽样。就算你推导出了条件分布的公式，你还需要有办法从这个分布中生成随机数。如果这些条件分布本身就是一些奇特的、没有标准生成算法的分布，那么你还是会卡住。

拒绝抽样 (Rejection Sampling)

核心思想：“接受/拒绝”机制

拒绝抽样的目标是从一个复杂的目标分布 $f(x)$ 中生成样本。它的方法不是直接生成，而是先从一个简单的候选分布 $g(x)$ 中“海量”地生成样本，然后按照特定规则进行筛选，只保留一部分样本，被保留下来的样本就符合我们想要的目标分布 $f(x)$。

• 目标分布 $f(x)$：我们想抽样的最终目标，但它很难直接处理。在图上是那条实线。
• 候选分布 $g(x)$：一个我们可以轻松抽样的简单分布，比如均匀分布或正态分布。
• 包络函数 $Cg(x)$：为了确保筛选的有效性，我们需要找到一个常数 $C$ ，使得 $C$ 倍的候选分布 $Cg(x)$ 能够完全“包住”目标分布 $f(x)$ 。

对候选分布 $g(x)$ 的要求

• $g(x)$ 必须易于抽样：这是整个方法的前提。如果我们连候选分布都很难抽样，那就毫无意义了。
• $Cg(x)$ 必须能“包住” $f(x)$ ：我们必须找到一个常数 $C$（通常 $C>1$ ），使得对于所有的 $x$ ，都有  $Cg(x) ≥ f(x)$ 。这个 $C$ 的最小值就是 $f(x)/g(x)$ 这个比值的上确界(supremum)。
• $Cg(x)$ 被称为包络函数 (envelope function)

![[Pasted image 20250908122247.png]]

算法步骤

• 概念版算法：
  1. 从候选分布 $g(x)$ 中抽取一个候选值 $x$。
  2. 以 $p=\frac{f(x)}{Cg(x)}$​ 的概率接受这个 $x$。
  3. 重复以上步骤直到获得足够多的样本。
• 实用版算法：
  1. 同时抽取两个随机数：一个候选值  $x\sim g(x)$  和一个均匀分布的随机数  $\mathrm{u\sim Unif}(0,1)$。
  2. 如果满足不等式 $u\leq\frac{f(x)}{Cg(x)}$，就接受这个 $x$。
  3. 重复以上步骤。

效率分析：$C$ 的含义与几何分布

几何分布说明了在一系列独立的伯努利试验中，为了首次获得成功，总共需要进行多少次试验。 $$P(Y=k)=(1-p)^{k-1}p,\quad k=1,2,\ldots,$$

• 如果每次试验的成功概率是p，那么在第k次试验才首次成功的概率是 P(Y=k) = (1-p)ᵏ⁻¹p。
• 它的期望（平均成功所需次数）是 1/p。

而在拒绝抽样算法中，我们为了获得一个被接受的目标样本，所需要生成的候选样本的总数量，就服从一个几何分布。

我们下面推导这个几何分布的“成功概率”，就是拒绝抽样的整体接受率。

一次抽样被“接受”的条件是 $u\leq\frac{f(x)}{Cg(x)}$

• 那么，对于一个固定的候选值 $x$，在 $u$ 很多的情况下， $u$ 的期望就是接受它的概率就是 $\frac{f(x)}{Cg(x)}$。
• 但是，候选值 $x$ 本身是随机的，它是从 $g(x)$ 中抽出来的。所以，$x$ 在不同位置出现的可能性不同。
• 为了得到总的、平均的成功概率，我们需要综合考虑所有可能的 $x$。在数学上，这就是求期望值 $E_g\left[\frac{f(x)}{Cg(x)}\right]$。
• 期望值的定义就是积分：∫(我们关心的量)×(该量出现的概率密度)dx。
• 所以，总成功率 = $\int\frac{f(x)}{Cg(x)}g(x)dx$。
• $g(x)$ 被约掉，积分变成了 $\int\frac{f(x)}{C}dx=\frac{1}{C}\int f(x)dx$。
• 因为 $f(x)$ 是一个概率密度函数，它在整个定义域上的积分必然等于1。所以总成功概率为 $1/C$。

那么，为了获得一个有效样本，平均需要生成的候选样本数量 = 几何分布的期望 = $1/p$ = $1/(1/C)$ = $C$。

这个结论告诉我们，包络常数 $C$ 不仅仅是一个抽象的数学符号，它有一个非常具体的物理意义：  $C$ 值，就代表了你平均要付出多少次计算（生成和检验候选样本），才能换来一个你真正想要的样本。

维度灾难

对于拒绝抽样，有以下五个要点：

• 关于包络常数 $C$ 的选择：任何比最优 $C$ 更大的 $C'$ 理论上都可行，但为了效率，我们应选择尽可能小的 $C$ 。
  • 任何比最优值C更大的常数C'在理论上都能让算法工作，但接受率是1/C'，所以C'越大，效率越低。我们应该让C'尽可能接近最优的C。
• 对数尺度计算：为了数值稳定性，实际编程中应在对数尺度下进行接受/拒绝的判断。
  • 所有操作都可以在对数尺度下进行：$log(u) ≤ log(f(x)) - log(g(x)) - log(C)$。
  • 这是因为概率值都是0到1之间的小数。当似然函数需要连乘上千个这样的数时，结果会迅速变得极小，超出计算机浮点数的表示精度，直接被当作0处理，导致计算错误。
  • 取对数 $log$ 可以将连乘变成连加，除法变成减法。l$og(a*b) = log(a)+log(b)$。这样，我们处理的都是一些绝对值比较正常的数字，计算过程会稳定得多。
• 维度灾难：拒绝抽样在高维空间中表现极差，此时必须使用MCMC等其他方法。
  • 在一维空间，我们是在一个面积里判断；在二维空间，我们是在一个体积里判断；在d维空间，我们是在一个d维的超体积里判断。
  • 随着维度的增加，“包络函数Cg(x)”和“目标函数f(x)”之间的“无效空间”会呈指数级增长。
  • 即使g(x)在每个维度上都和f(x)拟合得不错，但在高维空间中，一个点要同时在所有维度上都落在那个很小的目标区域内，其概率会变得微乎其微。
• 自适应 $C$：存在一些“自适应”的拒绝抽样变体，可以动态地优化 $C$。
  • 我们不预先固定一个Cg(x)包络，而是从一个粗糙的包络开始。在抽样过程中，我们利用那些被拒绝的点的信息，来不断地“修正”和“收紧”我们的包络线，让它越来越贴近f(x)，从而动态地提高抽样效率。
• 重尾条件：为了保证存在一个有限的 $C$，候选分布 $g(x)$ 的“尾巴”必须比目标分布 $f(x)$的“尾巴”更“重”。
  • $C$ 是否为无穷大，通常取决于候选分布的尾部，它必须比目标分布的尾部更“重”。
  • “重尾”的直观理解：当x趋向于无穷大时，概率密度函数f(x)下降到0的速度有多慢。下降得越慢，尾巴就越“重”。
  • 果候选分布g(x)的尾巴比f(x)“轻”（也就是g(x)下降到0的速度比f(x)快），那么在远离中心的地方，f(x)会“反超”g(x)。这样一来，比率f(x)/g(x)就会趋向于无穷大，我们就永远找不到一个有限的常数C来让Cg(x)能包住f(x)。为了确保安全，人们有时会选择一个已知尾巴很重的分布（如学生t分布）作为候选分布。

用先验直接抽样

目标分布是后验分布： $$\pi(\theta|y)\propto f(y|\theta)\pi(\theta)$$ 我们想从它采样。在拒绝采样里，我们需要一个“容易抽”的候选分布$g(\theta)$ ，一个最懒的想法是：直接把先验分布 $\mathrm\pi(\theta)$ 拿来当候选。那么就有$$\frac{\pi(\theta|y)}{C\pi(\theta)}=\frac{f(y|\theta)\pi(\theta)}{C\pi(\theta)}$$ $\mathrm\pi(\theta)$ 被约掉，接受概率简化成：$\frac{f(y|\theta)}{C}$ ，拒绝采样要求 $Cg(θ)$ 要“盖住”整个目标分布，即 $$C\geq f(y|\theta)\mathrm{~,~}\quad\forall\theta\mathrm{~}$$ 为了效率最高，我们取最小的 $C$，即：$$C^=\sup_\theta f(y|\theta)$$ 这就是似然函数的最大值。它在最大似然估计 (MLE) 的点 $\hat{\theta}$ 上取得。所以： $$C^=f(y|\hat{\theta})$$ 最终算法

• 从先验分布 $π(θ)$ 抽一个候选 $θ$。
• 从均匀分布 $u\sim\mathrm{Unif}[0,1]$ 抽一个随机数。
• 接受 $θ$ 如果：$$u\leq\frac{f(y|\theta)}{f(y|\hat{\theta})}$$

重要性抽样 (Importance Sampling)

核心思想：从“错误”的分布抽样，再用“权重”修正

当我们的目标分布 $f(x)$ 过于复杂，导致我们无法直接从中抽样时，重要性抽样提供了一套解决方案。

重要性权重：w = f/g

在贝叶斯推断中，我们最终得到了后验分布  $p(θ|D)$ 。但我们通常不只满足于知道这个分布的“形状”，我们更想计算它的一些统计特性，最常见的就是求期望 (Expectation)。

比如，我们想计算参数 $θ$ 的后验均值（即期望值）： $$E[\theta]=\int\theta\cdotp(\theta|D)d\theta$$ 或者计算某个关于 $θ$ 的函数 $h(θ)$ 的期望： $$E[h(\theta)]=\int h(\theta)\cdotp(\theta|D)d\theta$$ 这个积分通常非常复杂，无法直接用数学解析求解。一个很自然的蒙特卡洛想法是：

1. 从后验 $p(θ|D)$ 中抽取大量样本 $θ_1, θ_2, ..., θ_n$。
2. 然后用简单的平均来近似这个积分：$E[h(θ)] ≈ (1/n) * Σ h(θ_i)$。

这个想法的前提是你能轻易地从后验  $p(θ|D)$ 中抽样。但在现实中，后验分布往往形状奇特，直接从中抽样非常困难或低效（比如用拒绝采样，拒绝率可能极高）。

重要性抽样的思想是：

“既然从目标分布 $p(θ|D)$ 抽样那么难，我们能不能从一个我们熟悉的、简单的分布 $q(θ)$ 中抽样，然后通过某种方式，对抽出来的样本进行‘修正’，让它们看起来像是从 $p(θ|D)$ 中来的一样？”

这个“修正”的方式，就是为每个样本赋予一个“重要性权重” (Importance Weight)。

一个例子：民意调查

假设：

• 目标 p(θ|D)：你想知道全中国人民对某个问题的平均看法（比如，对甜粽子和咸粽子的喜爱程度）。这是你的目标分布。
• 困难：你没有资源去全国各地进行抽样调查。
• 简单分布 q(θ)：你只能在上海进行抽样调查（上海就是你的建议分布 q(θ)）。从上海市民中抽样很简单。

不能直接用上海的调查结果（样本均值）来代表全国的结果因为上海的样本分布 (q) 和全国的真实分布 (p) 存在偏差（比如，上海年轻人比例可能偏高，甜党可能占优等）。

重要性抽样的做法：

1. 在上海随机抽取了1000个市民（从q中抽样）。
2. 对每一个被抽到的市民，都计算一个权重 w。 w = 这个人在全国的代表性 / 这个人在上海的代表性 w = p(这个市民) / q(这个市民)
3. 如何计算权重？
   • 如果抽到了一个上海土生土长的年轻人。这种人在上海很常见（q值高），但在全国人口中的比例可能没那么高（p值相对低）。那么他的权重 w = p/q 就会小于1。我们需要降低他的发言权。
   • 如果抽到了一个来自农村的、在上海务工的老年人。这种人在上海可能比较稀有（q值低），但在全国人口中的比例其实不低（p值相对高）。那么他的权重 w = p/q 就会大于1。我们需要放大他的发言权。
4. 最终计算加权平均： 你用这1000个带权重的样本来计算全国的平均看法。 全国平均看法 ≈ Σ (每个人的看法 × 每个人的权重) / Σ (所有人的权重)

通过这个“加权修正”的过程，就巧妙地把一个有偏差的上海样本，变成了对全国情况的一个无偏估计。

比率形式 vs. 非比率形式

非比率形式

我们引入 $w_i=\frac{f(x_i)}{g(x_i)}$ 这个权重因子，那么对于我们的目标分布 $p(θ|D)$ ，我们需要知道其的归一化常数。假设我们知道，并且我们的 $p(θ|D)$ 是一个严格的、积分为1的概率密度函数（PDF）。

那么我们可以有以下的非比率形式 (Non-Ratio Form)： $$E[h(\theta)]\approx\frac{1}{N}\sum_{i=1}^Nw_i\cdot h(\theta_i)$$ 这里的权重 $w_i$ 是“真实权重”： $$w_i=\frac{p(\theta_i|D)}{q(\theta_i)}$$ 其中 $p(θ_i|D)$ 是真实的、已经归一化的后验概率密度。也就是说，$∫ p(θ|D) dθ = 1$。这个公式要求你率先知道后验分布的表达，这在实际上是不可能的，因为我们要求的就是后验分布的期望，我们本来就不知道后验分布是什么样的。 $$p(\theta|D)=\frac{p(\theta)\cdotp(D|\theta)}{p(D)}$$ 我们求的就是等号左边，它是未知的，因为 $p(D)$ 是个积分，我们求不出来。

非比率形式是一个理论上无偏（unbiased）的好估计量，但它要求我们知道一个我们实际上不知道的东西 ( $p(D)$ ），所以它在贝葉斯实践中几乎没有用武之地。

比率形式

这是我们在现实世界中真正使用的公式。我们面临的现实是：我们不知道 $p(θ|D)$，但我们能很轻松地获得它的未归一化形式，我们称之为 $p*(θ_i|D)$。

我们定义修正系数 = 这个点在真实分布中的重要性 / 这个点在抽样分布中的重要性

• 真实分布的重要性：由我们想模拟的目标 $p*(θ_i|D)$ 来衡量。
• 抽样分布的重要性：由我们实际使用的 $q(θ_i)$ 来衡量。 $$w_i^=\frac{p^(\theta_i|D)}{q(\theta_i)}$$ 这个时候我们对权重求加权平均 $$\text{加权平均}=\frac{\sum_iw_i\cdot v_i}{\sum_iw_i}$$ 我们想要估算的是 $h(θ)$ 的平均值。所以，对于我们抽到的每一个样本 $θ_i$ ，我们关心的“数值”就是 $h(θ_i)$。公式现在变成了： $$\text{期望}\approx\frac{\sum_iw_i\cdot h(\theta_i)}{\sum_iw_i}$$ 于是，我们得到了： $$E[h(\theta)]\approx\frac{\sum_iw_i^\cdot h(\theta_i)}{\sum_iw_i^}$$ 或者我们认为：
• 假设 $f(x)=C_f\cdot\tilde{f}(x)$ （我们只知道 $\tilde{f}(x)$）
• 假设 $g(x)=C_g\cdot\tilde{g}(x)$ （我们只知道 $\tilde{g}(x)$）
• 那么权重 $w_i=\frac{C_f\cdot\tilde{f}(x_i)}{C_g\cdot\tilde{g}(x_i)}$。这个权重里包含了一个我们不知道的常数$C_f/C_g$。

如果我们使用比率形式的估计： $$\frac{\sum w_ih(x_i)}{\sum w_i}=\frac{\sum\frac{C_f}{C_g}\tilde{w_i}h(x_i)}{\sum\frac{C_f}{C_g}\tilde{w_i}}$$ 常数约掉了，所以也是没问题的。所以在比率估计下期望值的估计是 $\frac{\sum w_ih(x_i)}{\sum w_i}$ ，这本质上是 $h(x_i)$ 的一个加权平均。

应用：方差缩减与敏感性分析

蒙特卡洛标准误差 (Monte Carlo standard error)

重要性抽样这个方法，最初被设计出来的目的，是为了提高蒙特卡洛积分的精度。当我们用蒙特卡洛方法估算一个值时（比如一个积分的期望），我们其实是在用一个随机样本的平均值来代替真实值。

如果运气不好，抽到了一批“差”的样本，估算结果可能就会离真实值很远。如果运气好，抽到一批“好”的样本，结果就可能很准。

每次重新做一次蒙特卡洛模拟，都会得到一个略微不同的答案。这种由于抽样的随机性而导致的估算结果的“抖动”或“不稳定性”，就用方差 (Variance) 或 标准误差 (Standard Error) 来衡量。方差越小，说明估算方法越稳定、越精确。用更少的样本就能得到更可靠的结果。

关于降低方差：

• 最简单粗暴的方法是增加样本量 n。样本越多，结果越稳定。但这会增加计算成本。
• 一个更聪明的方法，就是改进抽样策略本身。
• 重要性抽样的初衷，就是通过精心挑选建议分布 g(x)，来主动地、有策略地最小化最终估算结果的方差。

数学家已经证明，对于非比率形式的重要性抽样，能让最终估算结果方差最小的最优建议分布 $g_{optimal}(x)$，应该正比于： $$g_{optimal}(x)\propto|h(x)|\cdot f(x)$$

• f(x) 部分: 这很直观。我们的抽样重心 (g) 应该放在目标分布 f(x) 本身概率密度高的区域。即“在重要的地方多抽样”。
• |h(x)| 部分: 我们的抽样重心，还应该额外关注我们想求积分的那个函数 h(x) 绝对值大的区域。

注：

• 这个最优公式  $g ∝ |h(x)|f(x)$ 是针对非比率形式的重要性抽样推导出来的。
• 对于我们实际中更常用的比率形式，由于其结构更复杂（分母也是一个估计量），最小化其方差的最优  $g(x)$  的形式会不一样，也更复杂。

这个最优 $g(x)$ 的公式虽然完美，但在实践中却有一个悖论：

• 为了让 $g(x) = C |h(x)|f(x)$ 成为一个合法的概率分布，我们需要计算归一化常数 $C$ ，也就是计算积分 $\int|h(x)|f(x)dx$。
• 但这个积分往往和我们最初想解决的问题 $∫h(x)f(x)dx$ 一样困难，甚至更难。
• 虽然我们通常无法直接构造出这个“完美”的 $g(x)$ ，但这个公式为我们选择提议分布 $g(x)$ 提供了黄金指导原则：我们应该尽力选择一个我们已知且易于抽样的 $g(x)$ ，使其形状尽可能地模仿 $|h(x)|f(x)$ 的形状。

敏感性分析 (Sensitivity Analysis)

在贝叶斯分析中，选择先验分布 $π(θ)$ 是我们主观做出的决定。一个严谨的分析师必须回答一个问题：“我的最终结论，对我当初选择的这个先验有多敏感？”

• 如果我换一个稍微不同的先验，结论就天差地别，那说明我的结论非常不稳定，很大程度上是由我的主观臆断驱动的，而不是数据本身。
• 如果我换好几个合理的先验，结论都差不多，那说明我的结论是稳健的，主要是由数据驱动的。

敏感性分析就是做这个检查的过程。但如果每次检查都要重跑一次模型，成本太高了。

在这里，假如我们已经设定好了一个先验分布 $π(θ|ψ_0)$（其中 $ψ_0$ 是先验分布的超参数，比如Beta分布的 $α$ 和 $β$ ），并且通过一个非常耗时耗力的计算过程（比如MCMC模拟），已经从后验分布 $π(θ|y, ψ_0)$ 中获得了大量的样本 {$θ_i$}。

这时，我们或我们的合作者开始质疑：“如果当初我们的先验假设 $ψ_0$ 稍微不同，变成了 $ψ$ ，那结论会有什么变化呢？” 传统方法需要我们从头再跑一遍耗时耗力的模拟。

现在让我们回顾一下权重公式和贝叶斯定理： $$w_i=\frac{f(\theta_i)}{g(\theta_i)}=\frac{\pi(\theta_i|\mathbf{y},\psi)}{\pi(\theta_i|\mathbf{y},\psi_0)}$$ 根据贝叶斯定理，$后验 ∝ 似然 × 先验$，所以： $$\pi(\theta_i|\mathbf{y},\psi)\propto f(\mathbf{y}|\theta_i)\cdot\pi(\theta_i|\psi)$$ $$\pi(\theta_i|\mathbf{y},\psi_0)\propto f(\mathbf{y}|\theta_i)\cdot\pi(\theta_i|\psi_0)$$ 把它们代入到权重公式里： $$w_i\propto\frac{f(\mathbf{y}|\theta_i)\cdot\pi(\theta_i|\psi)}{f(\mathbf{y}|\theta_i)\cdot\pi(\theta_i|\psi_0)}$$ 代表着“数据证据”的似然函数 $f(y|θ_i)$ 在分子和分母中是完全一样的，所以可以直接约掉。最终的权重只取决于先验： $$w_i\propto\frac{\pi(\theta_i|\psi)}{\pi(\theta_i|\psi_0)}$$

覆盖

重要性抽样的前提是，建议分布 g（旧后验）必须在目标分布 f（新后验）所有重要的区域（即概率密度高的地方）都分配了不可忽略的概率。

假设旧后验 g 和新后验 f 都是山丘形状的曲线。

情况一：覆盖得很好 (Good Coverage)

• 旧后验 g（蓝色）和新后验 f（红色）的位置和形状都差不太多。
• 它们可能中心点稍微错开了一点，或者一个胖一点，一个瘦一点，但它们的主体部分是大量重叠的。
• 这意味着，f 这座红色山峰所在的位置，g 这座蓝色山丘也有相当的高度。

情况二：覆盖得极差 (Bad Coverage)

• 旧后验 g（蓝色）因为初始假设 ψ₀ 的原因，形成了一个在左边的山峰。
• 但你的新假设 ψ 导致新后验 f（红色）形成了一个在右边的山峰。
• 这两座山丘几乎没有任何重叠。红色山峰所在的位置，蓝色曲线的高度几乎是零。

因为重要性权重是 w(x) = f(x) / g(x)。

在上面那个“覆盖得极差”的图景中：

1. 绝大多数样本会怎样？
   • 我们的样本 {θᵢ} 是从 g（蓝色曲线）中抽取的，所以它们绝大多数都落在左边的蓝色山峰区域。
   • 在这些点上，g(θᵢ) 的值比较大，但 f(θᵢ) 的值几乎是零（因为红色山峰在右边）。
   • 所以，这些样本的权重 wᵢ = f(θᵢ)/g(θᵢ) = (几乎为0) / (比较大) ≈ 0。
   • 结果：大部分样本权重都接近于0，几乎成了“废样本”。
2. 极少数“幸运”样本会怎样？
   • g（蓝色曲线）可能会有一个非常非常长的尾巴，以至于它有极小的、非零的概率，抽到了一个正好落在右边红色山峰区域的样本，我们称之为 θ_lucky。
   • 在这个幸运点上，f(θ_lucky) 的值非常高（因为它在红色山峰的顶峰），而 g(θ_lucky) 的值极其微小（因为它在蓝色曲线遥远的尾部）。
   • 那么，这个幸运样本的权重 w_lucky = f(θ_lucky)/g(θ_lucky) = (非常大) / (极其小) = 很大的数
   • 结果：一个样本的权重可能会比其他所有样本的权重加起来还要大几百万倍。

例子

想象你是一名数据分析师，正在分析一个新网站按钮的点击率 θ。

1. 最初的分析 (旧模型):

• 初始信念 (旧先验)：根据过去的经验，这种按钮的点击率通常都在50%左右，不大可能太高或太低。所以，你选择了一个信息量很足的 (informative) 先验分布：Beta(α=10, β=10)。
  • 这个分布的形状是一个在 θ=0.5 处非常尖锐的山峰，它强烈地表达了“坚信点击率就在0.5附近”的观点。
• 结合了观测到的数据（比如1000次点击，有550次命中）和这个强大的先验，然后花费了大量计算资源（跑MCMC）得到了10000个后验样本 {θᵢ}。这些样本代表了在你的“专家经验”和“数据证据”结合后，你对点击率的最终认识。因为你的先验很强，这些样本点会紧密地聚集在0.5附近。

2. 新的疑问 (新模型):

• 老板说：“你的经验可能过时了。这次的按钮设计完全不同，我们应该假设我们对点击率一无所知，让数据自己说话。如果你当初用一个模糊的 (vague) 先验会怎样？”
• 老板建议的先验 (新先验)：Beta(α=1, β=1)。
  • 这个分布其实就是 Uniform(0, 1) 分布。它的形状是一条从0到1的平线，表达了“在看到数据前，我认为所有0到1之间的点击率都是完全等可能的”的开放心态。

在不重跑整个昂贵MCMC的前提下，快速回答老板的问题——如果当初用了Beta(1,1)这个“开放心态”的先验，我们对点击率的后验均值估计会是多少？

• 目标 f (新后验)：我们想知道但没有样本的后验分布 π(θ|y, ψ=(1,1))。
• 建议 g (旧后验)：我们已经有10000个样本的后验分布 π(θ|y, ψ₀=(10,10))。

只需要先计算权重 wᵢ

我们对已有的每一个样本 θᵢ，计算它的重要性权重。根据我们之前的推导，这个权重就是新旧先验的比值： $$w_{i} = \frac{新先验在 \theta_{i} 处的高度}{旧先验在 \theta_{i} 处的高度} = \frac{Beta(\theta_{i} | \alpha = 1, \beta = 1)}{Beta(\theta_{i} | \alpha = 10, \beta = 10)}$$

• 对于一个靠近0.5的样本 θᵢ (比如0.51):
  • 旧先验 Beta(10,10) 在这里的值非常高（因为它就是以0.5为中心的尖峰）。
  • 新先验 Beta(1,1) 在这里的值是1（因为它处处为1）。
  • 权重 wᵢ = 1 / (一个很高的值) = 一个很小的值。
  • 解读：这个样本点 θᵢ=0.51 在我们最初的模拟中被“过度”重视了，因为我们那个尖峰先验给了它特殊待遇。为了修正这一点，我们给它一个很低的权重，削弱它的影响力。
• 对于一个远离0.5的样本 θᵢ (比如0.8):
  • 旧先验 Beta(10,10) 在这里的值非常低（因为它在尖峰遥远的尾部）。
  • 新先验 Beta(1,1) 在这里的值还是1。
  • 权重 wᵢ = 1 / (一个很低的值) = 一个很大的值。
  • 解读：这个样本点 θᵢ=0.8 在我们最初的模拟中被“系统性地忽视”了。为了修正这一点，我们给它一个很高的权重，放大它的影响力，让它获得在“开放心态”下应有的话语权。

这个重新加权的过程，巧妙地“撤销”了旧先验（Beta(10,10)）所带来的偏见，并“应用”上了新先验（Beta(1,1)）的影响

最后，我们使用重要性抽样的比率公式来计算新的后验均值： $$\mathrm{New~Mean}\approx\frac{\sum_{i=1}^{10000}w_i\cdot\theta_i}{\sum_{i=1}^{10000}w_i}$$ 这个公式就是一个加权平均。我们把所有旧样本 {θᵢ} 拿过来，但不再把它们一视同仁地求平均，而是根据上面计算出的新权重 wᵢ 来求加权平均。那些在“开放心态”下更合理的样本，会在最终的均值计算中拥有更大的发言权。

抽样/重要性重抽样 (SIR)

拒绝抽样的缺点

拒绝抽样有两个主要痛点：

• 寻找一个最优（或足够好）的包络常数C是算法效率的关键，但这在实践中（尤其是在高维问题中）非常困难。
  • 拒绝抽样的生死命脉完全掌握在常数C的手中，因为接受率是1/C。
  • 为了让算法可行，我们必须保证C是f(x)/g(x)这个比值的上界。
  • 为了让算法高效，我们又希望C越小越好。
  • 在高维空间中，寻找这个“恰到好处”的C本身就是一个非常困难的优化问题，这使得拒绝抽样在很多复杂场景下难以应用。
• 拒绝抽样非常浪费，因为它会生成大量的“无用”样本，然后将它们直接丢弃。
  • 拒绝抽样的哲学是“非黑即白”：一个候选样本要么被100%接受，要么被100%丢弃。
  • 如果接受率只有5%，那意味着我们95%的计算（用于生成和评估那些被拒绝的样本）都被付诸东流了。

核心思想：结合重要性抽样和重抽样，避免浪费

面对上述困境，SIR方法提出了一种全新的、更“柔和”的思路。

SIR的核心思想是：我们从提议分布g(x)中生成的每一个样本，都或多或少地包含了关于目标分布f(x)的信息，不应该被粗暴地丢弃。

• 第一步：抽样/重要性加权 (Sampling/Importance)
  • 和重要性抽样一样，我们先从简单的提议分布 $g(x)$ 中生成一大批（比如 $N$ 个）候选样本。我们保留所有这些样本。
  • 然后，我们为每一个样本 $x_i$ 计算一个重要性权重 $w_i​=f(x_i​)/g(x_i​)$。这个权重衡量了该样本“有多么地像”一个来自目标分布 $f(x)$ 的样本。
• 第二步：重抽样 (Resampling)
  • 现在我们手上有一个带权重的样本集。我们可以把它看作是对目标分布 $f(x)$ 的一个离散近似。
  • 接下来，我们从这个加权的样本集中，进行有放回（或无放回）的加权随机抽样，来得到我们最终想要的、没有权重的 $m$ 个样本。
• 解决了 $C$ 的难题：在基础的SIR流程中，我们只需要计算每个点的权重比率 $f(x)/g(x)$ ，而不需要预先知道这个比率的全局最大值 $C$ 。这极大地降低了算法的使用门槛。
• 解决了“浪费”的问题：SIR没有丢弃任何一个初始样本。所有样本都被用来构建那个加权的离散分布，它们的信息通过权重得以保留。

算法步骤：“抽样-加权”与“重抽样”

SIR方法的目标是，从我们最终想研究的、但可能很复杂的目标分布 $f(x)$ 中，生成一个大小为 $m$ 的、近似独立同分布的样本集。

该方法分为两步：抽样步和重抽样步。

首先，我们从一个简单的提议分布 $g(x)$ 中抽取一个非常大的初始样本集，大小为 $J$（ $J ≥ m$ ）。然后，为这 $J$ 个样本中的每一个，都计算一个重要性权重 $w_j=f(X^{(j)})/g(X^{(j)})$ ，并将这些权重归一化，得到每个样本被选中的概率$ω_j$。

接着，我们从这 $J$ 个带权重的样本中，进行一次加权的、无放回的随机抽样，来抽取我们最终想要的 $m$ 个样本。

• SIR算法的效果，取决于两个因素：提议分布 $g(x)$ 与目标分布 $f(x)$ 的相似程度，以及初始样本量 $J$ 与最终样本量 $m$ 的比值 $J/m$ 的大小。
• $g(x)$ 越接近 $f(x)$ ，或者 $J/m$ 的比值越大，最终得到的样本质量就越好。

$J/m$ 比值的选择

我们第一步得到的那 $J$ 个带权重的样本，本质上是对连续的目标分布 $f(x)$ 的一个离散化的近似。

• 如果 $J$ 太小（比如 $J=m=100$ ），这个离散近似会非常粗糙。我们可能根本就没有抽到我们想要的样本，那么在第二步的重抽样中，需要的样本就永远不可能出现在我们的最终样本里，导致最终结果依然有偏。
• 如果 $J$ 很大（比如 $J=2000$ ），我们就有更大的机会在海选中覆盖到更多样化的人群（哪怕是那些在 $g(x)$ 中概率很低但在 $f(x)$ 中很重要的群体）。这样，我们构建的那个离散近似就更精确、更“高清”。在后续的重抽样中，我们就有了一个更优质的“候选池”来进行挑选。

理论上，当 $J/m$ 趋于无穷大时，SIR算法是精确的。实践中，文献建议 $J/m$ 的比值至少为10，比如20。

加权自助法 (Weighted Bootstrap)

抽样/重要性重抽样 (SIR) 有另一种、也是更常见的一种变体有放回的SIR (SIR with Replacement)。

有放回的SIR的关键区别就是在第二步的时候，要求我们从这个 $J$ 个样本池中，抽取 $m$ 个最终样本 ${Y^{(i)}}$ 。抽样方式需要满足一个条件：在期望（平均）上，第 $j$ 个候选样本 $X^{(j)}$ 在最终样本集里出现的次数 $q_j$ ，应该等于最终样本量 $m$ 乘以它的中奖概率 $ω_j$ ，即  $E(q_j) = mω_j$ 。

要满足上述抽象的数学条件，最简单直接的实现方式就是“简单的加权随机抽样（有放回）”，也就是所谓的加权自助法 (Weighted Bootstrap)。

在候选分布 $g$ 较差时，使用 带放回的重采样 更好。因为当当提议分布g和目标分布f相差很远时：

• 权重的极端化：重要性权重wⱼ = f(X⁽ʲ⁾)/g(X⁽ʲ⁾)会变得非常不稳定。绝大多数从g中抽到的样本，其f的值都很小，所以它们的权重wⱼ会趋近于0。而极少数幸运地落在f很高g很低的区域的样本，会得到极其巨大的权重。
• “黄金样本”的出现：我们的初始样本池J中，可能只有寥寥无几的几个“黄金样本”，而其他都是“炮灰样本”。

对比两种重抽样方式：

• 无放回 (Without Replacement)：在这种模式下，每个“黄金样本”最多只能被选中一次。如果我们的最终样本量m比“黄金样本”的数量还多，那么我们就不得不从那些权重很低的“炮灰样本”中进行凑数，这会污染我们最终样本的质量。
• 有放回 (With Replacement)：在这种模式下，一个“黄金样本”因为其极高的权重，可以被重复选中多次，这就像一个班里评选三好学生，某个同学实在太优秀了（权重极高），学校允许他一人独得三个名额。这使得我们的最终样本集m可以更多地由那些最能代表目标分布f的“黄金样本”构成，从而得到一个更准确的近似。