From 44e73b37e495045ce3024da93d4b74283f7a1f97 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sat, 18 Jul 2026 20:00:59 +1000 Subject: [PATCH 1/3] =?UTF-8?q?=F0=9F=94=84=20resync=20bayes=5Fnonconj.md?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .translate/state/bayes_nonconj.md.yml | 6 + lectures/bayes_nonconj.md | 912 +++++++++++--------------- 2 files changed, 402 insertions(+), 516 deletions(-) create mode 100644 .translate/state/bayes_nonconj.md.yml diff --git a/.translate/state/bayes_nonconj.md.yml b/.translate/state/bayes_nonconj.md.yml new file mode 100644 index 00000000..edbd635c --- /dev/null +++ b/.translate/state/bayes_nonconj.md.yml @@ -0,0 +1,6 @@ +source-sha: b78fbcddae98a645bff0f01bb28a1e7955db5f53 +synced-at: "2026-07-18" +model: claude-sonnet-5 +mode: RESYNC +section-count: 6 +tool-version: 0.17.0 diff --git a/lectures/bayes_nonconj.md b/lectures/bayes_nonconj.md index 14f1f286..d77c9f1f 100644 --- a/lectures/bayes_nonconj.md +++ b/lectures/bayes_nonconj.md @@ -9,655 +9,535 @@ kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 +translation: + title: 非共轭先验 + headings: + Overview: 概述 + The coin-flipping model: 抛硬币模型 + The coin-flipping model::Generating data: 生成数据 + The coin-flipping model::Specifying the model in NumPyro: 在NumPyro中设定模型 + MCMC reproduces the conjugate posterior: MCMC重现共轭后验 + Non-conjugate priors: 非共轭先验 + Non-conjugate priors::A uniform prior: 均匀先验 + Non-conjugate priors::A truncated log-normal prior: 截断对数正态先验 + Non-conjugate priors::A truncated Laplace prior: 截断拉普拉斯先验 + Variational inference: 变分推断 + Variational inference::Why variational inference?: 为什么要用变分推断？ + Variational inference::The evidence lower bound: 证据下界 + Variational inference::Implementing SVI in NumPyro: 在NumPyro中实现SVI + Variational inference::Comparing VI with MCMC: 将VI与MCMC进行比较 + Where to next: 接下来的方向 --- # 非共轭先验 -本讲是{doc}`quantecon讲座 `的续篇。 +```{include} _admonition/gpu.md +``` + +除了Anaconda中已有的库之外，本讲座还需要以下库： + +```{code-cell} ipython3 +:tags: [hide-output] -那节课在似然函数和参数先验分布恰好形成**共轭**对的情况下，提供了概率的贝叶斯解释，其中： +!pip install numpyro jax arviz +``` -- 应用贝叶斯法则产生的后验分布与先验具有相同的函数形式 +## 概述 -具有共轭关系的似然和先验可以简化后验的计算，有助于进行解析或近似解析计算。 +本讲是{doc}`prob_meaning`的续篇。 -但在许多情况下，似然和先验不需要形成共轭对。 +在那节课中，我们对硬币正面朝上的未知概率 $\theta$ 采用了**beta**先验分布，并配合**二项**似然函数。 -- 毕竟，一个人的先验是他或她自己的事情，只有在极小的巧合下才会采取与似然共轭的形式 -在这些情况下，计算后验概率会变得非常具有挑战性。 +该先验和似然构成了一对**共轭对**：应用贝叶斯法则得到的后验分布与先验属于*同一*分布族——同样是beta分布。 -在本讲中，我们将说明现代贝叶斯学者如何通过使用蒙特卡洛技术来处理非共轭先验，这涉及到： +共轭之所以便利，是因为它能给出闭式解的后验分布。 -- 首先巧妙地构建一个马尔可夫链，其不变分布就是我们想要的后验分布 -- 模拟该马尔可夫链直到其收敛，然后从不变分布中采样以近似后验分布 +但一个人的先验信念是他或她自己的事情，一般而言未必恰好与似然函数共轭。 -我们将通过使用两个强大的Python模块来说明这种方法，这些模块实现了这种方法以及下面将要描述的另一种密切相关的方法。 +当先验和似然**不**共轭时，后验通常没有闭式解，我们必须对其进行数值近似。 -这两个Python模块是： +本讲介绍两种被广泛使用的方法，二者都在概率编程库[NumPyro](https://num.pyro.ai/en/stable/getting_started.html)中实现： -- `numpyro` -- `pymc4` +* **马尔可夫链蒙特卡洛（MCMC）**——构造一个不变分布为后验分布的马尔可夫链，然后从中采样。我们使用**No-U-Turn采样器（NUTS）**，这是哈密顿蒙特卡洛的一种最先进形式。 -像往常一样，我们首先导入一些Python代码。 +* **变分推断（VI）**——用优化代替采样：在一族可处理的分布中搜索最接近后验分布的成员。 -```{code-cell} ipython3 -:tags: [hide-output] +(nuts)= +```{note} +本讲将NUTS视为一个黑箱。 -# install dependencies -!pip install numpyro pyro-ppl torch jax +简言之，它是**哈密顿蒙特卡洛**的一种形式，而哈密顿蒙特卡洛本身又是**Metropolis–Hastings**算法的一个版本：它提出候选抽样并接受或拒绝，从而使得到的马尔可夫链以后验分布为其不变分布。 + +它与基本的Metropolis–Hastings采样器的区别在于，其提议是根据对数后验的*梯度*（导数）信息构建的，这使得链能够高效地在参数空间中移动；此外NUTS还会自动调节每次提议移动的步长。 + +关于MCMC和Metropolis–Hastings算法更深入的介绍，参见[本讲](https://python-advanced.quantecon.org/mcmc.html)。 ``` + +我们的计划是： + +1. 确认MCMC能够重现我们可以解析计算的*共轭*beta后验分布——这可以在一个我们已知答案的问题上验证方法的有效性。 +2. 用若干**非共轭**先验代替beta先验，并用MCMC近似每个后验分布。 +3. 引入变分推断，并将其与MCMC进行比较。 + +让我们先导入一些库。 + ```{code-cell} ipython3 import numpy as np -import seaborn as sns import matplotlib.pyplot as plt -import matplotlib as mpl -FONTPATH = "fonts/SourceHanSerifSC-SemiBold.otf" -mpl.font_manager.fontManager.addfont(FONTPATH) -plt.rcParams['font.family'] = ['Source Han Serif SC'] - -from scipy.stats import binom import scipy.stats as st -import torch -# jax import jax.numpy as jnp -from jax import lax, random +from jax import random -# pyro -import pyro -from pyro import distributions as dist -import pyro.distributions.constraints as constraints -from pyro.infer import MCMC, NUTS, SVI, ELBO, Trace_ELBO -from pyro.optim import Adam - -# numpyro import numpyro -from numpyro import distributions as ndist -import numpyro.distributions.constraints as nconstraints -from numpyro.infer import MCMC as nMCMC -from numpyro.infer import NUTS as nNUTS -from numpyro.infer import SVI as nSVI -from numpyro.infer import ELBO as nELBO -from numpyro.infer import Trace_ELBO as nTrace_ELBO -from numpyro.optim import Adam as nAdam +import numpyro.distributions as dist +from numpyro.infer import MCMC, NUTS, SVI, Trace_ELBO +from numpyro.infer.autoguide import AutoNormal +from numpyro.optim import Adam + +import arviz as az ``` -## 在二项分布似然上释放MCMC -本讲座从{doc}`quantecon讲座`中的二项分布示例开始。 +## 抛硬币模型 -该讲座通过以下方式计算后验分布： +如{doc}`prob_meaning`中所述，硬币以概率 $\theta$ 正面朝上（$Y=1$），以概率 $1-\theta$ 反面朝上（$Y=0$）。 -- 通过选择共轭先验进行解析计算 +如果我们抛硬币 $n$ 次，正面朝上的次数 $k$ 服从**二项**分布 -本讲座则通过以下方式计算后验分布： +$$ +p(k \mid \theta) = \binom{n}{k}\, \theta^k (1-\theta)^{n-k} . +$$ -- 通过MCMC方法对后验分布进行数值采样，以及 -- 使用变分推断(VI)近似 +我们把 $\theta$ 视为具有先验密度 $p(\theta)$ 的随机变量，我们想要的是后验分布 -我们使用`pyro`和`numpyro`包，并借助`jax`来近似后验分布 +$$ +p(\theta \mid k) \propto p(k \mid \theta)\, p(\theta) . +$$ -我们使用几种不同的先验分布 +### 生成数据 -我们将计算得到的后验分布与{doc}`quantecon讲座`中描述的共轭先验相关的后验分布进行比较 +我们模拟一枚硬币的一系列抛掷结果，这枚硬币正面朝上的真实概率（分析者未知）为 $\theta = 0.4$。 +```{code-cell} ipython3 +def simulate_coin_flips(θ=0.4, n=20, seed=1234): + "抛硬币n次；返回由0（反面）和1（正面）组成的数组。" + rng = np.random.default_rng(seed) + return (rng.random(n) < θ).astype(int) + +data = simulate_coin_flips() +k, n = int(data.sum()), len(data) +k, n +``` -### 解析后验分布 +我们特意使用了一个**小**样本（$n = 20$）。 -假设随机变量$X\sim Binom\left(n,\theta\right)$。 +原因是先验分布的影响在数据稀少时最为显著。 -这定义了一个似然函数 +当样本量很大时，似然函数占主导地位，几乎任何合理的先验都会导致相同的后验分布——这正是我们在{doc}`prob_meaning`中看到的那种集中现象。 -$$ -L\left(Y\vert\theta\right) = \textrm{Prob}(X = k | \theta) = -\left(\frac{n!}{k! (n-k)!} \right) \theta^k (1-\theta)^{n-k} -$$ -其中 $Y=k$ 是一个观测数据点。 +适度的 $n$ 能保持先验的影响可见，这正是我们在此想要研究的。 -我们将 $\theta$ 视为一个随机变量，为其指定一个具有密度 $f(\theta)$ 的先验分布。 +### 在NumPyro中设定模型 -我们稍后会尝试其他先验分布，但现在，假设先验分布为 $\theta\sim Beta\left(\alpha,\beta\right)$，即： +对大多数读者来说，这将是第一次接触NumPyro，其风格需要一些时间来适应。 -$$ -f(\theta) = \textrm{Prob}(\theta) = \frac{\theta^{\alpha - 1} (1 - \theta)^{\beta - 1}}{B(\alpha, \beta)} -$$ +要使用它，我们把概率模型描述为一个Python函数——NumPyro有点令人困惑地把它称为**模型**（model）。 -我们现在选择这个作为先验分布，是因为我们知道二项分布似然函数的共轭先验是贝塔分布。 +这样的函数在被调用时不会*计算*任何东西，也不会返回后验分布。 -在 $N$ 个样本观测中观察到 $k$ 次成功后，$\theta$ 的后验概率分布为： +相反，它是对数据生成过程的一种*声明*：哪些量是随机的，它们服从什么分布，以及数据如何依赖于它们。 -$$ -\textrm{Prob}(\theta|k) = \frac{\textrm{Prob}(\theta,k)}{\textrm{Prob}(k)}=\frac{\textrm{Prob}(k|\theta)\textrm{Prob}(\theta)}{\textrm{Prob}(k)}=\frac{\textrm{Prob}(k|\theta) \textrm{Prob}(\theta)}{\int_0^1 \textrm{Prob}(k|\theta)\textrm{Prob}(\theta) d\theta} -$$ -=\frac{{N \choose k} (1 - \theta)^{N-k} \theta^k \frac{\theta^{\alpha - 1} (1 - \theta)^{\beta - 1}}{B(\alpha, \beta)}}{\int_0^1 {N \choose k} (1 - \theta)^{N-k} \theta^k\frac{\theta^{\alpha - 1} (1 - \theta)^{\beta - 1}}{B(\alpha, \beta)} d\theta} -$$ +一个推断算法——比如下面的NUTS采样器——随后会*读取*这个声明，并为我们求出后验分布。 -$$ -=\frac{(1 -\theta)^{\beta+N-k-1} \theta^{\alpha+k-1}}{\int_0^1 (1 - \theta)^{\beta+N-k-1} \theta^{\alpha+k-1} d\theta} . -$$ +在模型内部，每个随机量都通过调用`numpyro.sample`来引入，关键字`obs`决定了它的角色： + +* `numpyro.sample("θ", prior)` 引入一个名为`"θ"`的**潜在**（未观测）变量，从`prior`中抽取——这是我们希望推断的量。 + +* `numpyro.sample("k", dist.Binomial(n, θ), obs=k)` 引入一个**观测**变量：关键字`obs=k`将其固定为数据，这正是似然函数 $p(k \mid \theta)$ 进入模型的方式。 + +字符串名称（`"θ"`和`"k"`）是NumPyro用来跟踪这些变量的标签；我们稍后会用它们把后验抽样结果取出来。 + +我们只写*一个*模型，把先验分布作为参数传入，这样对于我们考虑的每一个先验——无论是否共轭——都可以原样复用它。 + +```{code-cell} ipython3 +def binomial_model(prior, k, n): + "带有调用者提供的θ先验的二项似然。" + θ = numpyro.sample("θ", prior) + numpyro.sample("k", dist.Binomial(n, θ), obs=k) +``` + +注意`binomial_model`不返回任何东西，而且我们从不自己调用它。 + +相反，我们把它交给一个推断算法，由算法提供参数并追踪这两条`sample`语句，从而组装出后验分布。 + +我们还编写一个小的辅助函数，用来在给定模型上运行NUTS并返回拟合好的采样器。 + +我们请求四条链，以便在下面检验收敛性，并用`chain_method="vectorized"`运行它们，这会在单个设备上同时评估所有链——因此同样的代码在CPU或GPU上都能不加修改地运行。 + +```{code-cell} ipython3 +def run_nuts(model, *args, seed=0, num_warmup=1000, num_samples=4000, num_chains=4): + "用NUTS采样器对一个NumPyro模型进行采样。" + mcmc = MCMC( + NUTS(model), + num_warmup=num_warmup, + num_samples=num_samples, + num_chains=num_chains, + chain_method="vectorized", + progress_bar=False, + ) + mcmc.run(random.key(seed), *args) + return mcmc +``` + +NumPyro建立在[JAX](https://docs.jax.dev)之上，而JAX显式地处理随机性：它不依赖全局随机状态，而是要求每次运行都拥有自己的**PRNG密钥**，这里通过`random.key(seed)`创建。 + +（这就是为什么我们在上面用NumPy的生成器来生成数据，而在这里使用JAX密钥的原因。） -因此， +`run_nuts`特意写得很通用：它对我们传入的任何模型进行采样，并通过`mcmc.run`把额外的参数（`*args`）转发给那个模型。我们始终以`run_nuts(binomial_model, prior, k, n)`的形式调用它，这样`prior`、`k`和`n`就会原样传递给`binomial_model`——自始至终只有这一个先验。 + +## MCMC重现共轭后验 + +在把MCMC用于困难问题之前，让我们先在一个简单问题上检验它。 + +对于$\text{Beta}(\alpha_0, \beta_0)$先验，后验分布可以解析求出（参见{doc}`prob_meaning`）： $$ -\textrm{Prob}(\theta|k) \sim {Beta}(\alpha + k, \beta+N-k) +\theta \mid k \sim \text{Beta}(\alpha_0 + k,\ \beta_0 + n - k) . $$ -以下Python代码实现了给定共轭beta先验的解析后验。 +我们取 $\alpha_0 = \beta_0 = 2$，并用NUTS对后验进行采样。 ```{code-cell} ipython3 -def simulate_draw(theta, n): - """ - 生成一个大小为n的伯努利样本，其中P(Y=1) = theta - """ - rand_draw = np.random.rand(n) - draw = (rand_draw < theta).astype(int) - return draw - - -def analytical_beta_posterior(data, alpha0, beta0): - """ - 给定观测数据，用参数(alpha, beta)的beta先验分布 - 解析计算后验分布 - - 参数 - --------- - num : int. - 计算后验时的观测数量 - alpha0, beta0 : float. - beta先验分布的参数 - - 返回值 - --------- - 后验beta分布 - """ - num = len(data) - up_num = data.sum() - down_num = num - up_num - return st.beta(alpha0 + up_num, beta0 + down_num) +α0, β0 = 2.0, 2.0 +mcmc = run_nuts(binomial_model, dist.Beta(α0, β0), k, n) ``` -### 近似后验分布的两种方法 -假设我们没有共轭先验。 +在查看后验之前，我们应当检查采样器是否已经完成了它的任务。 -那么我们就无法解析地计算后验分布。 +与我们习惯的独立抽样不同，MCMC返回的是一个*依赖*序列——一条马尔可夫链——其早期的抽样仍然记得链的起点。 -相反，我们使用计算工具来近似一组替代先验分布的后验分布，这需要用到Python中的`Pyro`和`Numpyro`包。 +只有当链已经"遗忘"其起点、并沉降到其不变分布——按构造正是我们想要的后验分布——之后，我们才能信任其输出。 -我们首先使用**马尔可夫链蒙特卡洛**（MCMC）算法。 +作为一种保障，我们从不同的随机起点运行了**四条**链（`run_nuts`中的`num_chains=4`），现在检验它们是否彼此一致。 -我们实现NUTS采样器来从后验分布中采样。 +[ArviZ](https://www.arviz.org/)是一个用于检验贝叶斯采样器输出的配套库。 -通过这种方式，我们构建一个近似后验分布的采样分布。 +函数`az.from_numpyro`把我们的NumPyro结果重新打包为ArviZ的标准数据结构，`az.summary`则打印出每个参数的汇总统计和收敛性诊断表。 -在此之后，我们部署另一个称为**变分推断**（VI）的程序。 +```{code-cell} ipython3 +idata = az.from_numpyro(mcmc) +az.summary(idata, var_names=["θ"]) +``` -特别是，我们在`Pyro`和`Numpyro`中都实现了随机变分推断（SVI）机制。 +这张表中有两列是值得理解的收敛性诊断指标。 -MCMC算法据说能产生更准确的近似，因为原则上它直接从后验分布中采样。 +* **`r_hat`**（Gelman–Rubin统计量）比较每条链*内部*的抽样离散程度与各条链*之间*的离散程度。如果所有链都已经收敛到同一分布，这两者应当相符，`r_hat`接近$1.0$；若数值超过大约$1.01$，则说明各链彼此不一致，抽样结果尚不可信。 -但是它在计算上可能很昂贵，尤其是当维度很大时。 -VI方法可能更便宜，但很可能会产生较差的后验近似，原因很简单，因为它需要猜测一个用于近似后验的参数化**指导函数形式**。 +* **`ess_bulk`**和**`ess_tail`**报告*有效样本量*。由于连续的MCMC抽样是相关的，一条长度为$N$的链所携带的信息少于$N$个独立抽样所能携带的信息；有效样本量估计的是它相当于多少个独立抽样（分别在分布的主体部分和尾部）。数值越大越好。 -这个指导函数充其量也只能是一个不完美的近似。 +这里`r_hat`基本上是$1.0$，有效样本量达到数千，说明各条链的混合情况良好。 -通过限制假定后验具有受限函数形式所付出的代价，后验近似问题被转化为一个明确的优化问题，该问题寻求假定后验的参数，以最小化真实后验和假定后验分布之间的Kullback-Leibler (KL)散度。 +**迹图（trace plot）**给出了同一事实的直观检验。 - - 最小化KL散度等价于最大化一个称为**证据下界**（ELBO）的标准，我们很快就会验证这一点。 +ArviZ的`plot_trace`为每个参数绘制两个面板：右侧是抽样值随迭代次数的变化（每条链一种颜色的线）；左侧是各条链抽样的密度估计。 -## 先验分布 +混合良好的链在右侧看起来像平稳噪声——一条模糊而扁平的带状区域，各条链彼此重叠而非漂移或游走——它们在左侧的密度曲线几乎完全重合。 -为了能够应用MCMC采样或VI，`Pyro`和`Numpyro`要求先验分布满足特殊性质： -- 我们必须能够从中进行采样； -- 我们必须能够逐点计算对数概率密度函数； -- 概率密度函数必须对参数可微。 +```{code-cell} ipython3 +az.plot_trace(idata, var_names=["θ"]) +plt.tight_layout() +plt.show() +``` -我们需要定义一个分布`class`。 +我们的链通过了这两项检验，因此可以信任这些抽样结果，转而来看后验分布本身。 -我们将使用以下先验： +现在我们把MCMC后验与解析beta后验进行比较。 -- 在区间$[\underline \theta, \overline \theta]$上的均匀分布，其中$0 \leq \underline \theta < \overline \theta \leq 1$。 +```{code-cell} ipython3 +θ_grid = np.linspace(0.001, 0.999, 500) +samples = np.asarray(mcmc.get_samples()["θ"]) + +fig, ax = plt.subplots() +ax.hist(samples, bins=50, density=True, alpha=0.4, + label="MCMC后验") +ax.plot(θ_grid, st.beta(α0 + k, β0 + n - k).pdf(θ_grid), + 'k-', lw=2, label="解析后验") +ax.plot(θ_grid, st.beta(α0, β0).pdf(θ_grid), + 'C1--', lw=2, label="先验") +ax.set_xlabel(r"$\theta$") +ax.legend() +plt.show() +``` -- 支撑在$[0,1]$上的截断对数正态分布，参数为$(\mu,\sigma)$。 +MCMC抽样的直方图恰好落在解析后验密度之上。 - - 要实现这一点，令$Z\sim Normal(\mu,\sigma)$且$\tilde{Z}$为支撑在$[\log(0),\log(1)]$上的截断正态分布，则$\exp(Z)$具有支撑在$[0,1]$上的对数正态分布。这很容易编码，因为`Numpyro`内置了截断正态分布，而`Torch`提供了包含指数变换的`TransformedDistribution`类。 -- 另外，我们可以使用拒绝采样策略，将界限外的概率率设为$0$，并通过原始分布的CDF计算的总概率来重新缩放被接受的样本（即在界限内的实现值）。这可以通过使用`pyro`的`dist.Rejector`类来定义截断分布类来实现。 +采样器行之有效，因此我们可以在没有闭式解后验的先验上依赖它。 - - 我们在下面的部分实现这两种方法，并验证它们产生相同的结果。 +## 非共轭先验 -- 一个支撑限制在$[0,1]$区间内的偏移冯·米塞斯分布，其参数为$(\mu,\kappa)$。 +我们现在保持二项似然和同样的数据不变，但把beta先验替换为与之**不**共轭的先验。 - - 设$X\sim vonMises(0,\kappa)$。我们知道$X$的支撑范围是$[-\pi, \pi]$。我们可以定义一个偏移的冯·米塞斯随机变量$\tilde{X}=a+bX$，其中$a=0.5, b=1/(2 \pi)$，这样$\tilde{X}$的支撑范围就在$[0,1]$上。 +对每一个先验，配方都是一样的： - - 这可以使用`Torch`的`TransformedDistribution`类及其`AffineTransform`方法来实现。 -- 如果我们想要先验服从冯·米塞斯分布(von-Mises)且中心为$\mu=0.5$,我们可以选择一个较高的集中度参数$\kappa$,使得大部分概率质量位于$0$和$1$之间。然后我们可以使用上述策略进行截断。这可以通过`pyro`的`dist.Rejector`类来实现。在这种情况下,我们选择$\kappa > 40$。 +1. 描述该先验，并将其构建为一个NumPyro分布， +2. 把它传入`binomial_model`并运行NUTS， +3. 把先验和得到的后验画在一起。 -- 一个截断的拉普拉斯分布。 +下面这个辅助函数在同一坐标轴上画出先验密度和后验抽样。 - - 我们还考虑了截断的拉普拉斯分布,因为它的密度函数呈现分段非光滑的形式,并具有独特的尖峰形状。 +```{code-cell} ipython3 +def plot_prior_posterior(prior, samples, title=""): + "在[0, 1]上叠加绘制θ的先验密度和后验MCMC抽样。" + grid = jnp.linspace(0.001, 0.999, 500) + # 将密度限制在先验的支撑范围内：dist.Uniform.log_prob + # 即使在[low, high]之外也会返回其常数值 + in_support = np.asarray(prior.support(grid)) + prior_pdf = np.where(in_support, np.exp(np.asarray(prior.log_prob(grid))), 0.0) + + fig, ax = plt.subplots() + ax.hist(np.asarray(samples), bins=50, density=True, alpha=0.4, + label="后验（MCMC）") + ax.plot(np.asarray(grid), prior_pdf, 'C1--', lw=2, label="先验") + ax.set_xlabel(r"$\theta$") + ax.set_xlim(0, 1) + ax.legend() + if title: + ax.set_title(title) + plt.show() +``` - - 可以使用`Numpyro`的`TruncatedDistribution`类创建截断的拉普拉斯分布。 +### 均匀先验 + +最简单的非共轭先验是**均匀**先验：分析者认为某个区间内$\theta$的每个取值都同等可能。 + +在整个$[0, 1]$区间上取均匀先验表示无差异态度。 + +由于其密度为常数，此时后验分布仅与似然函数成正比。 ```{code-cell} ipython3 -# 由Numpyro使用 -def TruncatedLogNormal_trans(loc, scale): - """ - 使用numpyro的TruncatedNormal和ExpTransform获取截断对数正态分布 - """ - base_dist = ndist.TruncatedNormal(low=jnp.log(0), high=jnp.log(1), loc=loc, scale=scale) - return ndist.TransformedDistribution( - base_dist,ndist.transforms.ExpTransform() - ) - -def ShiftedVonMises(kappa): - """ - 使用AffineTransform获取平移的冯·米塞斯分布 - """ - base_dist = ndist.VonMises(0, kappa) - return ndist.TransformedDistribution( - base_dist, ndist.transforms.AffineTransform(loc=0.5, scale=1/(2*jnp.pi)) - ) - -def TruncatedLaplace(loc, scale): - """ - 获取区间[0,1]上的截断拉普拉斯分布 - """ - base_dist = ndist.Laplace(loc, scale) - return ndist.TruncatedDistribution( - base_dist, low=0.0, high=1.0 - ) +mcmc_flat = run_nuts(binomial_model, dist.Uniform(0.0, 1.0), k, n) +plot_prior_posterior(dist.Uniform(0.0, 1.0), + mcmc_flat.get_samples()["θ"], + title="平坦均匀先验") +``` + +后验分布集中在样本频率 $k/n$ 附近，正如似然函数所显示的那样。 + +现在假设分析者反而确信这枚硬币偏向正面，于是在$[0.5, 0.95]$上取均匀先验。 + +这个先验对真实值$\theta = 0.4$附近的区域赋予了*零*密度。 + +```{code-cell} ipython3 +mcmc_restr = run_nuts(binomial_model, dist.Uniform(0.5, 0.95), k, n) +plot_prior_posterior(dist.Uniform(0.5, 0.95), + mcmc_restr.get_samples()["θ"], + title="限制性均匀先验") +``` + +后验分布无法在先验为零的地方分配概率质量，因此它堆积在下边界$0.5$附近——尽可能靠近先验所允许的数据方向。 + +这是一个生动的警示：一个排除了真相的先验，无论收集多少数据都永远无法被推翻。 + +### 截断对数正态先验 + +均匀先验是平坦的。更现实的先验则是光滑且不对称的。 -# 由Pyro使用 -class TruncatedLogNormal(dist.Rejector): - """ - 通过Pyro中的拒绝采样定义截断对数正态分布 - """ - def __init__(self, loc, scale_0, upp=1): - self.upp = upp - propose = dist.LogNormal(loc, scale_0) - - def log_prob_accept(x): - return (x < upp).type_as(x).log() - - log_scale = dist.LogNormal(loc, scale_0).cdf(torch.as_tensor(upp)).log() - super(TruncatedLogNormal, self).__init__(propose, log_prob_accept, log_scale) - - @constraints.dependent_property - def support(self): - return constraints.interval(0, self.upp) - - -class TruncatedvonMises(dist.Rejector): - """ - 通过Pyro中的拒绝采样定义截断冯·米塞斯分布 - """ - def __init__(self, kappa, mu=0.5, low=0.0, upp=1.0): - self.low, self.upp = low, upp - propose = dist.VonMises(mu, kappa) - - def log_prob_accept(x): - return ((x > low) & (x < upp)).type_as(x).log() - - log_scale = torch.log( - torch.tensor( - st.vonmises(kappa=kappa, loc=mu).cdf(upp) - - st.vonmises(kappa=kappa, loc=mu).cdf(low)) - ) - super(TruncatedvonMises, self).__init__(propose, log_prob_accept, log_scale) - - @constraints.dependent_property - def support(self): - return constraints.interval(self.low, self.upp) +在$[0, 1]$上一个方便的选择是**截断对数正态**分布：取$Z \sim N(\mu, \sigma)$并截断到$Z \le 0$，令$\theta = e^{Z}$，这样它就落在$(0, 1]$内。 + +NumPyro通过让`TruncatedNormal`经过`ExpTransform`来构造这个分布。 + +```{code-cell} ipython3 +def truncated_lognormal(μ, σ): + "截断到单位区间(0, 1]的对数正态分布。" + base = dist.TruncatedNormal(loc=μ, scale=σ, low=-jnp.inf, high=0.0) + return dist.TransformedDistribution(base, dist.transforms.ExpTransform()) + +prior_ln = truncated_lognormal(0.0, 1.0) +mcmc_ln = run_nuts(binomial_model, prior_ln, k, n) +plot_prior_posterior(prior_ln, mcmc_ln.get_samples()["θ"], + title="截断对数正态先验") ``` -### 变分推断 -变分推断方法不直接从后验分布中采样，而是用一族可处理的分布/密度来近似未知的后验分布。 +该先验偏好较小的$\theta$值，但由于$\sigma = 1$使其相当发散，似然函数把后验拉向样本频率。 + +我们保留`mcmc_ln`——下面将把它与变分推断进行比较。 + +### 截断拉普拉斯先验 + +我们最后一个先验有一个尖锐的、非光滑的峰值。 + +**拉普拉斯**密度 $\propto e^{-|\theta - \mu| / b}$ 在其中心$\mu$处有一个拐点，表示一种强烈的信念：$\theta$位于$\mu$附近，同时仍允许尾部出现意外。 + +我们把它截断到$[0, 1]$上，并以$0.5$为中心。 + +```{code-cell} ipython3 +def truncated_laplace(μ, b): + "截断到单位区间[0, 1]的拉普拉斯分布。" + return dist.TruncatedDistribution(dist.Laplace(μ, b), low=0.0, high=1.0) + +prior_lp = truncated_laplace(0.5, 0.1) +mcmc_lp = run_nuts(binomial_model, prior_lp, k, n) +plot_prior_posterior(prior_lp, mcmc_lp.get_samples()["θ"], + title="截断拉普拉斯先验") +``` + +这个带尖峰的先验把后验拉向$0.5$，偏离了接近$0.4$的样本频率。 + +这里的拉力比较温和，因为先验虽然有峰值，但并不十分陡峭；如果$b$更小，它就会主导这个规模不大的样本。 + +NUTS无需任何特殊调节即可处理先验中的这个拐点——这是基于梯度的采样器搭配自动微分的一个实际优势。 + +## 变分推断 -然后，它寻求最小化近似分布与真实后验分布之间的统计差异度量。 +MCMC通过从后验中*采样*来近似后验。 -因此，变分推断(VI)通过求解最小化问题来近似后验分布。 +**变分推断（VI）**采取了不同的路径：它把后验近似转化为一个*优化*问题。 -设我们要推断的潜在参数/变量为$\theta$。 +我们把注意力限制在一族可处理的密度 $q_\phi(\theta)$——**引导分布（guide）**——上，它由参数$\phi$索引，我们搜索该族中最接近后验的成员。 -设先验分布为$p(\theta)$，似然函数为$p\left(Y\vert\theta\right)$。 +### 为什么要用变分推断？ -我们想要求得$p\left(\theta\vert Y\right)$。 +如果NUTS已经能返回精确的后验，为什么还要引入另一种方法？ -根据贝叶斯法则： +答案是**规模**。 + +MCMC在每一步都要在整个数据集上评估似然函数，而且它所需要的步数往往随参数的维度增长。 + +对于大型数据集或高维模型——例如机器学习中常见的层级模型和神经网络——这可能会慢到不切实际的程度。 + +变分推断的伸缩性要好得多，因为其目标函数（下面将介绍的ELBO）可以用在小型随机数据子集上计算的*随机*梯度来最大化——这与训练深度学习模型所用的机制相同。 + +它还能给出一个紧凑的参数化近似，事后存储和抽样的成本都很低。 + +代价是精度：VI只能返回*引导分布族内*的最佳拟合，并且可能低估不确定性。 + +一个经验法则是：当你需要精确的后验且问题规模小到可以承受时，优先选择MCMC；当模型对MCMC来说太大，或者一个快速的近似答案已经足够好时，选择VI。 + +### 证据下界 + +设先验为$p(\theta)$，似然为$p(Y \mid \theta)$，其中$Y$表示观测数据（这里是正面次数$k$）。 + +根据贝叶斯法则， $$ -p\left(\theta\vert Y\right)=\frac{p\left(Y,\theta\right)}{p\left(Y\right)}=\frac{p\left(Y\vert\theta\right)p\left(\theta\right)}{p\left(Y\right)} +p(\theta \mid Y) = \frac{p(Y, \theta)}{p(Y)} = \frac{p(Y \mid \theta)\, p(\theta)}{p(Y)}, $$ 其中 $$ -p\left(Y\right)=\int d\theta p\left(Y\mid\theta\right)p\left(Y\right). +p(Y) = \int p(Y \mid \theta)\, p(\theta)\, d\theta . $$ (eq:intchallenge) -{eq}`eq:intchallenge`右侧的积分通常很难计算。 -考虑一个由参数$\phi$参数化的**引导分布**$q_{\phi}(\theta)$，我们将用它来近似后验分布。 +{eq}`eq:intchallenge`中的积分是麻烦所在：在非共轭情形下它没有闭式解。 -我们选择引导分布的参数$\phi$，以最小化近似后验分布$q_{\phi}(\theta)$与后验分布之间的Kullback-Leibler (KL)散度： +我们用**Kullback–Leibler（KL）散度**来度量引导分布$q_\phi(\theta)$与后验分布之间的差异 $$ - D_{KL}(q(\theta;\phi)\;\|\;p(\theta\mid Y)) \equiv -\int d\theta q(\theta;\phi)\log\frac{p(\theta\mid Y)}{q(\theta;\phi)} +D_{KL}\big(q_\phi(\theta)\ \|\ p(\theta \mid Y)\big) += -\int q_\phi(\theta)\, \log \frac{p(\theta \mid Y)}{q_\phi(\theta)}\, d\theta , $$ -因此，我们需要一个能解决以下问题的**变分分布**$q$： +并选择$\phi$使其最小化。 + +KL散度仍然涉及难以处理的后验分布，但我们可以对其进行重新整理。利用$p(\theta \mid Y) = p(\theta, Y) / p(Y)$， $$ -\min_{\phi}\quad D_{KL}(q(\theta;\phi)\;\|\;p(\theta\mid Y)) +\begin{aligned} +D_{KL}\big(q_\phi \,\|\, p(\theta \mid Y)\big) + & = -\int q_\phi(\theta)\, \log \frac{p(\theta, Y) / p(Y)}{q_\phi(\theta)}\, d\theta \\ + & = -\int q_\phi(\theta) \left[\log \frac{p(\theta, Y)}{q_\phi(\theta)} - \log p(Y)\right] d\theta \\ + & = -\int q_\phi(\theta)\, \log \frac{p(\theta, Y)}{q_\phi(\theta)}\, d\theta + \log p(Y) , +\end{aligned} $$ -注意到： +其中最后一行用到了$\int q_\phi(\theta)\, d\theta = 1$。整理后得到， $$ -\begin{aligned}D_{KL}(q(\theta;\phi)\;\|\;p(\theta\mid Y)) & =-\int d\theta q(\theta;\phi)\log\frac{P(\theta\mid Y)}{q(\theta;\phi)}\\ - & =-\int d\theta q(\theta)\log\frac{\frac{p(\theta,Y)}{p(Y)}}{q(\theta)}\\ - & =-\int d\theta q(\theta)\log\frac{p(\theta,Y)}{p(\theta)q(Y)}\\ - & =-\int d\theta q(\theta)\left[\log\frac{p(\theta,Y)}{q(\theta)}-\log p(Y)\right]\\ -$$ -& =-\int d\theta q(\theta)\log\frac{p(\theta,Y)}{q(\theta)}+\int d\theta q(\theta)\log p(Y)\\ - & =-\int d\theta q(\theta)\log\frac{p(\theta,Y)}{q(\theta)}+\log p(Y)\\ -\log p(Y)&=D_{KL}(q(\theta;\phi)\;\|\;p(\theta\mid Y))+\int d\theta q_{\phi}(\theta)\log\frac{p(\theta,Y)}{q_{\phi}(\theta)} -\end{aligned} +\log p(Y) = D_{KL}\big(q_\phi \,\|\, p(\theta \mid Y)\big) + + \underbrace{\int q_\phi(\theta)\, \log \frac{p(\theta, Y)}{q_\phi(\theta)}\, d\theta}_{\text{ELBO}} . $$ -对于观测数据$Y$，$p(\theta,Y)$是一个常数，所以最小化KL散度等价于最大化 +左边的边际似然$\log p(Y)$不依赖于$\phi$。 + +因此**最小化**KL散度等价于**最大化**第二项，即**证据下界（ELBO）**： $$ -ELBO\equiv\int d\theta q_{\phi}(\theta)\log\frac{p(\theta,Y)}{q_{\phi}(\theta)}=\mathbb{E}_{q_{\phi}(\theta)}\left[\log p(\theta,Y)-\log q_{\phi}(\theta)\right] +\text{ELBO}(\phi) \equiv \int q_\phi(\theta)\, \log \frac{p(\theta, Y)}{q_\phi(\theta)}\, d\theta += \mathbb{E}_{q_\phi(\theta)}\big[\log p(\theta, Y) - \log q_\phi(\theta)\big] . $$ (eq:ELBO) -公式{eq}`eq:ELBO`被称为证据下界(ELBO)。 - -可以使用标准优化程序来搜索我们参数化分布$q_{\phi}(\theta)$中的最优$\phi$。 +由于$D_{KL} \ge 0$，ELBO是$\log p(Y)$的一个下界——这也是其名称的由来。 -参数化分布$q_{\phi}(\theta)$被称为**变分分布**。 -我们可以在Pyro和Numpyro中使用`Adam`梯度下降算法来实现随机变分推断(SVI)以近似后验分布。 +关键在于，{eq}`eq:ELBO`只涉及*联合*密度 $p(\theta, Y) = p(Y \mid \theta)\, p(\theta)$，这是我们可以计算的，而不涉及难以处理的归一化常数$p(Y)$。 -我们使用两组变分分布：Beta分布和支撑在$[0,1]$上的截断正态分布 +这个期望可以通过从$q_\phi$中采样来估计，而$\phi$可以通过梯度上升来改进——这就是**随机变分推断（SVI）**。 - - Beta分布的可学习参数是(alpha, beta)，两者都是正数。 - - 截断正态分布的可学习参数是(loc, scale)。 +### 在NumPyro中实现SVI -我们将截断正态分布的'loc'参数限制在区间$[0,1]$内。 +我们需要一个引导分布$q_\phi$。 -## 实现 +最简单的选择是**自动引导（autoguide）**：NumPyro检查模型并自动为我们构建一个引导分布。 -我们构建了一个Python类`BaysianInference`，初始化时需要以下参数： +`AutoNormal`在每个潜在变量上放置一个独立的正态分布，并经过变换以满足其支撑范围——这里是为了把$\theta$保持在$(0, 1)$内。 -- `param`：依赖于分布类型的参数元组/标量 -- `name_dist`：指定分布名称的字符串 +我们把SVI应用到上面的截断对数正态模型上，并用Adam优化器最大化ELBO。 -(`param`, `name_dist`)配对包括： -- ('beta', alpha, beta) - -- ('uniform', upper_bound, lower_bound) +```{code-cell} ipython3 +guide = AutoNormal(binomial_model) +optimizer = Adam(step_size=0.01) +svi = SVI(binomial_model, guide, optimizer, loss=Trace_ELBO()) -- ('lognormal', loc, scale) - - 注意：这是截断的对数正态分布。 -- ('vonMises', kappa)，其中kappa表示集中参数，中心位置设为$0.5$。 - - 注意：在使用`Pyro`时，这是原始vonMises分布的截断版本； - - 注意：在使用`Numpyro`时，这是**平移后**的分布。 +svi_result = svi.run(random.key(0), 5000, prior_ln, k, n, progress_bar=False) +``` -- ('laplace', loc, scale) - - 注意：这是截断的拉普拉斯分布 +SVI最大化ELBO；等价地，它最小化其负值，即所报告的损失。 -类`BaysianInference`有几个关键方法： -- `sample_prior`: - - 可用于从给定的先验分布中抽取单个样本。 +一条趋于平坦的损失曲线表明已经收敛。 -- `show_prior`: - - 通过重复抽样并拟合核密度曲线来绘制近似的先验分布。 +```{code-cell} ipython3 +fig, ax = plt.subplots() +ax.plot(svi_result.losses) +ax.set_xlabel("步数") +ax.set_ylabel("负ELBO") +ax.set_title("SVI收敛情况") +plt.show() +``` -- `MCMC_sampling`: - - 输入：(data, num_samples, num_warmup=1000) - - 接收一个`np.array`数据并生成大小为`num_samples`的后验MCMC采样。 +### 将VI与MCMC进行比较 -- `SVI_run`: - - 输入：(data, guide_dist, n_steps=10000) - - guide_dist = 'normal' - 使用**截断的**正态分布作为参数化的guide -- guide_dist = 'beta' - 使用beta分布作为参数化的指导分布 - - 返回值: (params, losses) - 以`dict`形式存储的学习参数和每一步的损失向量。 +为了评估这一近似，我们从拟合好的引导分布中抽样，并将其与同一（对数正态先验）模型的NUTS后验进行比较。 ```{code-cell} ipython3 -class BayesianInference: - def __init__(self, param, name_dist, solver): - """ - 参数 - --------- - param : tuple. - 包含分布所有相关参数的元组对象 - dist : str. - 分布的名称 - 'beta', 'uniform', 'lognormal', 'vonMises', 'tent' - solver : str. - pyro或numpyro - """ - self.param = param - self.name_dist = name_dist - self.solver = solver - - # jax需要显式传入PRNG状态 - self.rng_key = random.PRNGKey(0) - - - def sample_prior(self): - """ - 定义在Pyro/Numpyro模型中用于采样的先验分布。 - """ - if self.name_dist=='beta': - # 解包参数 - alpha0, beta0 = self.param - if self.solver=='pyro': - sample = pyro.sample('theta', dist.Beta(alpha0, beta0)) - else: - sample = numpyro.sample('theta', ndist.Beta(alpha0, beta0), rng_key=self.rng_key) - - elif self.name_dist=='uniform': - # 解包参数 - lb, ub = self.param - if self.solver=='pyro': - sample = pyro.sample('theta', dist.Uniform(lb, ub)) - else: - sample = numpyro.sample('theta', ndist.Uniform(lb, ub), rng_key=self.rng_key) - - elif self.name_dist=='lognormal': - # 解包参数 - loc, scale = self.param - if self.solver=='pyro': - sample = pyro.sample('theta', TruncatedLogNormal(loc, scale)) - else: - sample = numpyro.sample('theta', TruncatedLogNormal_trans(loc, scale), rng_key=self.rng_key) - - elif self.name_dist=='vonMises': - # 解包参数 - kappa = self.param - if self.solver=='pyro': - sample = pyro.sample('theta', TruncatedvonMises(kappa)) - else: - sample = numpyro.sample('theta', ShiftedVonMises(kappa), rng_key=self.rng_key) - - elif self.name_dist=='laplace': - # 解包参数 - loc, scale = self.param - if self.solver=='pyro': - print("警告：请使用Numpyro进行截断拉普拉斯分布。") - sample = None - else: - sample = numpyro.sample('theta', TruncatedLaplace(loc, scale), rng_key=self.rng_key) - - return sample - - - def show_prior(self, size=1e5, bins=20, disp_plot=1): - """ - 通过从先验分布采样并绘制近似采样分布来可视化先验分布 - """ - self.bins = bins - - if self.solver=='pyro': - with pyro.plate('show_prior', size=size): - sample = self.sample_prior() - # 转换为numpy - sample_array = sample.numpy() - - elif self.solver=='numpyro': - with numpyro.plate('show_prior', size=size): - sample = self.sample_prior() - # 转换为numpy - sample_array=jnp.asarray(sample) - - # 绘制直方图和核密度估计 - if disp_plot==1: - sns.displot(sample_array, kde=True, stat='density', bins=bins, height=5, aspect=1.5) - plt.xlim(0, 1) - plt.show() - else: - return sample_array - - - def model(self, data): - """ - 通过指定先验分布、条件似然和数据条件来定义概率模型 - """ - if not torch.is_tensor(data): - data = torch.tensor(data) - # 设置先验 - theta = self.sample_prior() - - # 从条件似然中采样 - if self.solver=='pyro': - output = pyro.sample('obs', dist.Binomial(len(data), theta), obs=torch.sum(data)) - else: - # 注意：numpyro.sample()要求obs=np.ndarray - output = numpyro.sample('obs', ndist.Binomial(len(data), theta), obs=torch.sum(data).numpy()) - return output - - - def MCMC_sampling(self, data, num_samples, num_warmup=1000): - """ - 使用MCMC数值计算给定数据下的后验分布，先验为由(alpha0, beta0)参数化的beta分布 - """ - # 使用pyro - if self.solver=='pyro': - # 张量化 - data = torch.tensor(data) - nuts_kernel = NUTS(self.model) - mcmc = MCMC(nuts_kernel, num_samples=num_samples, warmup_steps=num_warmup, disable_progbar=True) - mcmc.run(data) - - # 使用numpyro - elif self.solver=='numpyro': - data = np.array(data, dtype=float) - nuts_kernel = nNUTS(self.model) - mcmc = nMCMC(nuts_kernel, num_samples=num_samples, num_warmup=num_warmup, progress_bar=False) - mcmc.run(self.rng_key, data=data) - - # 收集样本 - samples = mcmc.get_samples()['theta'] - return samples - - - def beta_guide(self, data): - """ - 定义用于在Pyro/Numpyro中近似后验的候选参数化变分分布 - 这里我们使用参数化beta分布 - """ - if self.solver=='pyro': - alpha_q = pyro.param('alpha_q', torch.tensor(0.5), - constraint=constraints.positive) - beta_q = pyro.param('beta_q', torch.tensor(0.5), - constraint=constraints.positive) - pyro.sample('theta', dist.Beta(alpha_q, beta_q)) - - else: - alpha_q = numpyro.param('alpha_q', 10, - constraint=nconstraints.positive) - beta_q = numpyro.param('beta_q', 10, - constraint=nconstraints.positive) - - numpyro.sample('theta', ndist.Beta(alpha_q, beta_q)) - - - def truncnormal_guide(self, data): - """ - 定义用于在Pyro/Numpyro中近似后验的候选参数化变分分布 - 这里我们使用[0,1]上的截断正态分布 - """ - loc = numpyro.param('loc', 0.5, - constraint=nconstraints.interval(0.0, 1.0)) - scale = numpyro.param('scale', 1, - constraint=nconstraints.positive) - numpyro.sample('theta', ndist.TruncatedNormal(loc, scale, low=0.0, high=1.0)) - - - def SVI_init(self, guide_dist, lr=0.0005): - """ - 使用Adam优化器初始化SVI训练模式 - 注意：truncnormal_guide只能与numpyro求解器一起使用 - """ - adam_params = {"lr": lr} - - if guide_dist=='beta': - if self.solver=='pyro': - optimizer = Adam(adam_params) - svi = SVI(self.model, self.beta_guide, optimizer, loss=Trace_ELBO()) - - elif self.solver=='numpyro': - optimizer = nAdam(step_size=lr) - svi = nSVI(self.model, self.beta_guide, optimizer, loss=nTrace_ELBO()) - - elif guide_dist=='normal': - # 仅允许numpyro - if self.solver=='pyro': - print("警告：请使用Numpyro和TruncatedNormal指导") - svi = None - - elif self.solver=='numpyro': - optimizer = nAdam(step_size=lr) - svi = nSVI(self.model, self.truncnormal_guide, optimizer, loss=nTrace_ELBO()) - else: - print("警告：请输入'beta'或'normal'") - svi = None - - return svi - - def SVI_run(self, data, guide_dist, n_steps=10000): - """ - 运行SVI并返回优化后的参数和损失 - - 返回值 - -------- - params : 指导分布的学习参数 - losses : 每一步的损失向量 - """ - - # 初始化SVI - svi = self.SVI_init(guide_dist=guide_dist) - - # 执行梯度步骤 - if self.solver=='pyro': - # 张量化数据 - if not torch.is_tensor(data): - data = torch.tensor(data) - # 存储损失向量 - losses = np.zeros(n_steps) - for step in range(n_steps): - losses[step] = svi.step(data) - - # pyro仅支持beta VI分布 - params = { - 'alpha_q': pyro.param('alpha_q').item(), - 'beta_q': pyro.param('beta_q').item() - } - - elif self.solver=='numpyro': - data = np.array(data, dtype=float) - result = svi.run(self.rng_key, n_steps, data, progress_bar=False) - params = dict( - (key, np.asarray(value)) for key, value in result.params.items() - ) - losses = np.asarray(result.losses) - - return params, losses +vi_samples = guide.sample_posterior( + random.key(1), svi_result.params, sample_shape=(4000,) +)["θ"] +nuts_samples = mcmc_ln.get_samples()["θ"] + +fig, ax = plt.subplots() +ax.hist(np.asarray(nuts_samples), bins=50, density=True, alpha=0.4, + label="MCMC（NUTS）") +ax.hist(np.asarray(vi_samples), bins=50, density=True, alpha=0.4, + label="VI（AutoNormal）") +ax.set_xlabel(r"$\theta$") +ax.legend() +plt.show() ``` + +两种近似在后验分布的位置和离散程度上大致一致。 + +它们不必完全一致。 + +MCMC对真实后验进行采样（存在蒙特卡洛误差），而VI给出的是*其引导分布族内*的最佳拟合。 + +平均场正态引导分布在变换后的尺度上是对称的，可能会漏掉真实后验中的偏度或重尾。 + +这是成本与精度之间的权衡：VI用优化代替采样，在高维情形下往往快得多，但其近似质量的上限取决于引导分布的灵活性。 + +## 接下来的方向 + +本讲展示了当先验和似然不共轭时，如何利用NumPyro中的NUTS和随机变分推断来计算后验分布。 + +同样的工具也可以用于更丰富的模型。 + +{doc}`ar1_bayes`和{doc}`ar1_turningpts`两讲将NumPyro应用于自回归时间序列的贝叶斯估计和预测，其中参数是一个向量，无法进行共轭分析。 \ No newline at end of file From 4c19212fc6aedb40211c54237affd8069da18fcd Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sat, 18 Jul 2026 21:20:44 +1000 Subject: [PATCH 2/3] Restore CJK font configuration dropped by the resync The resynced file keeps Chinese plot labels but reverted the import cell to the source's exact form, losing the Source Han Serif setup - Chinese in figures would render as missing glyphs. Residual #107 class found in the merge review; applied wave-wide. Co-Authored-By: Claude Fable 5 --- lectures/bayes_nonconj.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lectures/bayes_nonconj.md b/lectures/bayes_nonconj.md index d77c9f1f..5b9b952a 100644 --- a/lectures/bayes_nonconj.md +++ b/lectures/bayes_nonconj.md @@ -84,6 +84,10 @@ translation: ```{code-cell} ipython3 import numpy as np import matplotlib.pyplot as plt +import matplotlib as mpl +FONTPATH = "fonts/SourceHanSerifSC-SemiBold.otf" +mpl.font_manager.fontManager.addfont(FONTPATH) +plt.rcParams['font.family'] = ['Source Han Serif SC'] import scipy.stats as st import jax.numpy as jnp From db492a62843a585ebc0126eaf7166b00eb453188 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sun, 19 Jul 2026 10:32:46 +1000 Subject: [PATCH 3/3] Restore trailing newline (engine issue tracked in QuantEcon/action-translation#116) Co-Authored-By: Claude Fable 5 --- lectures/bayes_nonconj.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/bayes_nonconj.md b/lectures/bayes_nonconj.md index 5b9b952a..14d98360 100644 --- a/lectures/bayes_nonconj.md +++ b/lectures/bayes_nonconj.md @@ -544,4 +544,4 @@ MCMC对真实后验进行采样（存在蒙特卡洛误差），而VI给出的 同样的工具也可以用于更丰富的模型。 -{doc}`ar1_bayes`和{doc}`ar1_turningpts`两讲将NumPyro应用于自回归时间序列的贝叶斯估计和预测，其中参数是一个向量，无法进行共轭分析。 \ No newline at end of file +{doc}`ar1_bayes`和{doc}`ar1_turningpts`两讲将NumPyro应用于自回归时间序列的贝叶斯估计和预测，其中参数是一个向量，无法进行共轭分析。