贝叶斯数据分析（4）——单参数模型（三）

本笔记基于Bayesian Data Analysis Third Edition， Andrew Gelman et. al. 学习编写，由于是英文教材，可能在学习过程中一些翻译或者内容有误，如有问题或错误，可以发送至邮箱[email protected]反馈。在学习本书中，需要有一定数学基础或者数理统计的基础，并且存在很多需要计算的场景，对于每一条定理，从头开始的证明会让你更明白每一步是如何实现的。

那么，我们讲完了很多模型的设置，一般为了方便，我们会选用共轭先验来进行计算。在这种情况下，我们需要知道的是，我们假设先验的时候，其可以包括很多信息，但同时，它也可以几乎没有信息（因为我们无法观察），也可以介乎两者之间。但我们之前说过，后验是先验和数据的折衷，数据少时，先验会明显影响后验；数据多时，似然通常会主导后验。

我们这里选用一个例子讲述有信息的先验：

有信息先验（Informative prior distribution）

背景：美国各县在 1980-1989 年间白人男性肾脏/输尿管癌死亡率。原始地图显示：死亡率最高的 10% 县很多在美国中部；但再看死亡率最低的 10% 县，也有很多在美国中部。

如果我们熟悉美国地理的就知道，美国的常住人口基本位于东部和西部，而中部的居民比较少（由于环境不如东西部好），因此我们真的可以说中部地区XX县的肾癌死亡率高或者低吗，显然不可以，因为中部有很多人口很少的县，且肾癌本身处于罕见疾病，小县在十年内可能出现 0 例死亡，也可能偶然出现 1 例死亡。由于分母很小，1 例死亡就会让原始死亡率显得很高；0 例死亡又会让原始死亡率显得很低。接下来我们构建模型来做一示范。需要注意的是，由于这些地区之前本身有异质性，因此使用分层（hierarchal model ）会更加合适，但为了简单，我们暂时不考虑。

Poisson-Gamma 模型

对于罕见疾病（比如癌症），我们通畅会使用泊松分布来进行拟合发病率。

对每个县 $j$，设：

• $y_j$：1980-1989 年间该县观察到的肾癌死亡人数；
• $n_j$：该县人口数；
• $\omega_j$：该县年死亡率，单位是 deaths per person-year；
• $10n_j$：十年人年数，也就是 person-years。

于是我们便可以写出似然：

而对于总体的肾癌死亡率为$\frac{y_j}{10n_j}$

为了估计每个县的死亡率，我们使用Gamma先验

在这个参数下，根据Gamma分布的性质，我们写下期望和方差：

在这里，我们可以通过设置不同的 $\alpha, \beta$ 来进行先验的构造，在这里，我们需要明白的是，我们构造的先验是没有唯一正确的，在这里我们使用 $\alpha=20,\quad \beta=430000$ 来构造先验。

此时：

在这里，这个期望由县级死亡率数据的总体特征来近似确定。换句话说，我们把原有的信息纳入了我们的模型。

后验分布

此时，我们就可以计算后验分布：

此时，后验均值和方差为：

我们令暴露量：$E_j=10n_j$， 则一般形式下，我们可以计算一个普遍式的后验均值：

我们可以改写：

我们可以发现：这个后验均值就是先验和数据的折衷：

后验均值 = 先验权重 * 先验均值 + 数据权重 * 死亡率

换句话说，当我们人口越小的县，那么其 $E_j$ 越小，数据权重越低，后验估计越接近先验均值。人口越多的县，$E_j$ 越大，数据权重越高，后验估计越接近该县的原始死亡率。

先验预测分布

$y_j$ 的边际分布为：

我们之前也说过，此时计算的结果是负二项分布：

其中，均值和方差为：

我们看回方差，我们可以知道其第一项来自 Poisson 抽样波动，第二项来自县与县之间死亡率 $\omega_j$ 的差异。

预测

我们首先看真实情况：

我们可以清晰的看到人口小的县更容易出现极端死亡率。

小县人口

我们假设一个人口数少的县，n=1,000

那么，十年的人年数为10,000:

如果，在这十年内没有人死亡，即$y_j=0$，死亡率是 $0$，后验均值为：

当 $y_j=1$，死亡率为 1/ 10000，即 $10^{-4}$，此时，后验均值为：

当 $y_j=2$，死亡率是 2/10000， 此时后验均值为：

对小县来说，0 例、1 例、2 例造成的死亡率差异非常大，但贝叶斯估计会收缩到总体先验均值附近。这种收缩是在承认小样本条件下原始率噪声很大的基础上，把总体信息纳入估计。

大县人口

如果 $n_j=1000000$，则 $E_j=10n_j=10000000$. 在进行先验预测模拟后，我们得到：十年死亡人数的中位数约为 473，50% 区间约为 $[393,545]$

假设，当 $y_j=393$ 时，其死亡率为 $3.93\times 10^{-5}$ ，则后验均值为：

而当 $y_j=545$ 时，死亡率为 $5.45\times 10^{-5}$ , 后验均值为：

对于大县来说，真实死亡率和后验均值非常接近，因为数据量已经足够大，先验的影响很小

下面这幅图就很明显展示不同人口下其死亡率的估计，且随着人口数增大，其置信区间也会变小。

构造信息先验

我们先计算可观测的死亡率：

此时对于先验预测，我们计算其均值和方差为：

由于不同县的 $n_j$ 不同，实际估计时可以用样本中 $1/(10n_j)$ 的平均值来近似上式中的暴露量项，再把观测到的死亡率均值、方差与理论均值、方差匹配，解出 $\alpha,\beta$。这是一种经验贝叶斯或矩匹配思路：用总体的数据帮助构造先验，同时是对分层模型的近似。

非信息先验（Non-informative prior distribution）

合适的先验与不合适的先验

其实是Proper prior和improper prior，更多的是它是否在后验收敛，为了减少犯错概率后面用英文代称。

如果一个先验密度 $p(\theta)$ 不依赖数据，并且满足：

那么我们称作为proper prior。其实通常的，只要其是收敛于一个常数即可，当不收敛于1而收敛于其他常数时，我们称作非正态密度，但是其可以被标准化，通过乘以积分值倒数即可。而如果其不收敛：

则它是 improper prior。improper prior 不是合法的概率分布。但需要注意的是，我们仍然可以用improper prior得到proper posterior。因此二者并不是等价的。

我们看回已知方差的正态分布的例子：

我们之前算过，此时后验的precision为：

我们知道，当先验的precision远小于数据本身的precision时，我们会有：

此时，我们等价于将先验里的 $\tau_0^2 \to \infin$ ，则此时，$p(\theta)\propto 1,\qquad -\infty<\theta<\infty$

那么显然，$\int p(\theta)\,d\theta=\infty$，其是一个improper prior，但注意到，只要我们有一个观测值 $n > 0$ 即有proper posterior。

因此，既然improper prior可以算出proper posterior，那么prior的设定是否proper好像就不重要了，但事实情况是，我们要明白我们的目的是获得未知参数 $\theta$ 的分布，因此我们想要的是后验分布，所以后验是否proper是更加重要的，因此，我们必须检查的是：

如果这个积分有限，就可以把未归一化后验归一化为 proper posterior；如果积分无限，后验本身就不是概率分布，那么后续的推断就没有意义。

Jeffreys 不变性原则和先验

不变性原则

在选择非信息先验的时候，我们首先要知道我们选择的分布，参数都是人为选择的，这种情况下会引入主观性，因此Jeffreys提出了解决方案：首先我们在这里先借用一个分布，即二项分布，以二项分布中的成功概率为例，参数可以直接写成：$\theta\in(0,1)$ 同时也可以写成 logit 尺度：

此时，$\theta$ 和 $\phi$ 表示同一个未知成功概率，但所在的尺度不同。非信息先验希望尽量减少人为尺度选择对推断的影响。在这种要求下，使用 $\theta$ 作为参数和使用 $\phi=h(\theta)$ 作为参数时，同一套先验构造规则应当给出彼此等价的先验密度。

我们用一个例子来看：

在均匀先验下：$p_\theta(\theta)=1,\qquad 0<\theta<1$

若令：$\phi=\operatorname{logit}(\theta)$

则：

根据密度变换公式（一维Jacobian）:

因此：

这个密度在 $\phi$ 上是随 $\phi$ 变化的函数。概率尺度上的均匀先验变换到 logit 尺度后，会变成一个中间高、两端低的密度。

假设我们在 logit 尺度上取均匀先验：$p_\phi(\phi)\propto 1$

则变换回概率尺度有

这对应 improper 的 $\operatorname{Beta}(0,0)$ 形式。因此，我们发现 $\theta$ 尺度上的均匀先验 $\operatorname{Beta}(1,1)$ 与 logit 尺度上的均匀先验二者不等价。

Jeffreys prior

给定似然：$p(y\mid \theta)$, 记对数似然为：$\ell(\theta)=\log p(y\mid \theta)$

参数 $\theta$ 的 Fisher 信息定义为：

则Jeffreys 先验定义为;

$J(\theta)$ 衡量数据对参数 $\theta$ 的局部辨识能力。Fisher 信息越大，似然在该参数值附近对参数变化越敏感（二阶导）；Fisher 信息越小，似然在该区域越平缓。Jeffreys 先验使用 $\sqrt{J(\theta)}$ 作为密度形状，使先验构造规则在一一参数变换下保持一致。但应该注意的是，其仍然是一个先验设定，其特殊之处在于，先验密度由抽样模型的 Fisher 信息推出，避免选定一个主观分布族再指定超参数。某些模型中的 Jeffreys 先验是 proper prior，某些模型中的 Jeffreys 先验仍然是 improper prior，仍需检查后验密度是否收敛。

书上的主要内容如上，但我觉得对于我而言这是一个比较难以理解的知识点

Jeffreys prior 服从 Jeffreys不变性原则

接下来我们进行一个证明：

设参数的一一变换为 $\phi=h(\theta),\quad\ \theta=h^{-1}(\phi)$.

我们对似然进行对数转化，则有：

由链式法则：

两边平方并取期望：

即：

另一方面，密度变换公式为：

若在 $\theta$ 尺度上使用 Jeffreys 先验：

当变换到 $\phi$ 尺度下时：

二项分布的Jeffreys prior distribution

依旧我们使用二项分布来为我们的例子：

则似然为：

取对数：

计算对数似然的二阶导数：

于是乎，计算Fisher信息：

因此Jeffreys 先验为：

由于 $n$ 与 $\theta$ 无关，可以并入归一化常数：

所以我们构造的先验为：

所以后验为：

仍然为Beta分布，$\theta\mid y\sim \operatorname{Beta}\left(y+\frac12,n-y+\frac12\right)$， 根据Beta分布性质，当$\alpha>0, \quad \beta>0$， 这个后验积分是有限的，因此为proper posterior。

枢轴量（Pivotal quantities）

定义:设总体 $X$ 有概率密度（或分布律）$f(x;\theta)$ ，其中 $\theta$ 是待估的未知参数。设 $X_1,\dots,X_n$ 是一个样本，记：

$$G = G(X_1,\dots,X_n;\theta)$$

为样本和待估参数 $\theta$ 的函数，如果 $G$ 的分布已知，不依赖与任何参数，就称 $G$ 为枢轴量

也就是说，枢轴量与某个待估参数有关（其主要用于未知参数的区间估计），且本身含有未知参数（待估参数），因此具有“不可观察性”，也就是说即使选定了样本也无法计算出确定的值；其分布是明确的（有具体的数学公式，不包含未知参数）

例如：若有 $y\mid \theta\sim N(\theta,\sigma^2)$ , 其中 $\sigma^2$ 已知，且设 $u=y-\theta$ ，在给定 $\theta$ (即条件化)， $u\mid \theta=y-\theta\sim N(0,\sigma^2)$ ，此时$u$ 的分布不含 $\theta$，因此 $u=y-\theta$ 是一个枢轴量。

枢轴量意义在于：把数据相对于参数的偏离标准化成一个分布已知的量。在构造非信息先验时，可以要求在观察到数据后，这个枢轴量的后验分布仍然保持其原来的抽样分布。这个要求的成立假设为：在缺少额外先验信息时，数据只改变参数的位置或尺度，不额外改变枢轴量本身的分布形状。

对于位置参数，参数只控制分布的平移。并且在没有额外信息时，一个自然的做法是在位置尺度上不给任何位置更多权重，因此得到：

对于尺度参数，参数控制分布的伸缩。在没有额外信息时，一个自然做法是对乘法尺度上的相对变化给相同权重，也就是在 $\log\theta$ 上取均匀先验，因此得到：

枢轴量在构造非信息先验的适用范围有限。它适合纯位置参数和纯尺度参数这类结构清楚的单参数模型。对于多参数模型、分层模型或参数之间相互耦合的模型，枢轴量可能不唯一，依据不同枢轴量得到的先验也可能不同。

非信息先验局限

1. 如果数据真的足够强，那么在一组相对平坦的先验之间，不论哪一个选择通常不会显著影响后验
2. 平坦性依赖参数化。一个在 $\theta$ 上均匀的先验，在 $\log\theta$ 或 $\theta^2$ 上就不再均匀
3. improper prior 在模型比较或模型平均中会造成额外困难，因为不同模型中的归一化常数可能没有定义

弱信息先验 （Weak-informative prior distribution）

弱信息先验是 proper prior，其放入了少量的信息，它的目的是提供足够约束，使后验不会跑到明显不合理的区域。因为在很多实际问题中，我们并不是对先验分布一无所知，即便不能精确写出完整的先验分布，也通常知道参数不可能离谱到某些程度。

也就是说，弱信息先验是介于强信息先验和非信息先验的居中，我们举个例子：

对于性别而言，我们不知道这真正的比例，但我们知道这个比例不太可能远离0.5， 因此此时可以使用集中在 0.4 到 0.6 附近的先验（或者更加严格），例如：$N(0.5,0.1^2)$，如果直接把正态分布用于概率参数，可以把它理解为边界外概率极小的近似，或在 $[0,1]$ 上截断。为了保持二项模型的共轭形式，也可以使用：$\theta\sim \operatorname{Beta}(20,20)$，80% 的概率会落在0.4-0.6. 这个先验很松，但其排除了明显不合理的范围。

在构造弱信息先验时，我们主要可以从两点出发：

1. 从非信息先验出发，通过加上一点先验信息，使得推断在合理范围以内。比如估计罕见病患病率，我们不能使用一个 $\theta\sim U(0,1)$， 且其真实患病率可能小于 1/10000，因为假设我们在100个人里观测到0个人，那么此时后验为 $\theta\mid y\sim \operatorname{Beta}(1,101)$, 则此时后验均值为 1/102，显然远大于 1/10000. 所以不太符合
2. 从强信息先验出发，把它适当放宽，以反映历史信息和当前研究场景之间可能存在的差异。比如历史数据可能强烈提示某个治疗效应为正，但当前实验面对的人群、剂量、测量方式可能不同。我们可以保留效应大致为正的信息，但扩大先验方差，避免历史信息过度支配新数据。

在有些场景下，即使我们想要的确实是方向性先验，也可能不希望把它直接放进分析中。比如假设我们相信某个处理效应为正，但实验本身是为了检验这个理论。如果直接使用偏向正效应的先验，就会降低证据门槛。在这种情况下，我们可以使用关于 0 对称的弱信息先验，如 $\theta\sim N(0,A^2)$ ，这样既能限制极端不合理的效应大小，又不会预先把结论推向正方向或负方向。

总结