如果移动端无法正确渲染某些公式,建议转PC站点,也是做了自适应处理的。点击这里跳转

BMSI系列的第一篇,介绍了贝叶斯框架下的单参数模型,同时也算是贝叶斯入门的第一步。

手敲$\LaTeX$难免出现纰漏,有任何疑似错误或者不清楚的地方请直接在下方评论区留言(评论区只开在【Series】下),谢谢各位读者。

一个简单的引入——二项分布

模型假设

我们从二项分布中抽取了一系列样本$x_1,\ldots,x_n$,我们想知道这个二项分布的参数$\theta$。

传统统计学

可以轻松得出,无论是极大似然估计还是矩估计都会给出结果$\theta=y/n$,其中$y$为$x_1,\ldots,x_n$中为1的个数。

这个结果是相当自然而符合直觉的,但其也有我们不满意的地方。一个浅显的例子就是当我们进行2次实验得到2次1,和我们进行20000次实验得到20000次1,我们得到的参数估计都是$1$.但直觉上后者为1的可信度就远远高于前者。区间估计可以刻画出这样的估计的一个“可信度”,但对于点估计,我们只能给出$1$这个结论。

与之相比,贝叶斯框架下给出的点估计就看起来更加合理,我们在下面进行简单的介绍。

贝叶斯方法

正如贝叶斯学派所认为的,模型的参数并不是“确定但隐藏”的,而是服从一个分布——我们称之为先验(prior)$p(\theta)$。先验就“先”在我们此时还没有任何的数据,只能根据我们的常识或者信念给出这样的分布。与此同时,我们对模型有所假设,因此可以得到概率密度函数$p(y\lvert \theta)$,这刻画着当前的参数对这组数据的解释度。使用贝叶斯公式

$$p(\theta,y)=p(\theta)p(y\lvert \theta)=p(y)p(\theta\lvert y)$$

我们可以得到参数的后验分布(posterior)

$$p(\theta\lvert y)=p(\theta)p(y\lvert \theta)/p(y)\propto p(\theta)p(y\lvert \theta)$$

具体到当前的例子。由于我们对参数的分布没有任何先验的认识,我们采取$\left[0,1\right]$上的均匀分布作为先验——也就是我们认为$\theta$此时的取值都是等可能的,我们并不对任何取值有偏向性。(PS:与之相对的,对于抛硬币的问题,我们可能会更倾向于参数为$1/2$,因此我们选择的先验可能为一个单峰且峰在$1/2$处的一个分布)

而概率密度函数则就是简单的二项分布,即

$$p(y\lvert \theta)=\binom{n}{y}\theta^y (1-\theta)^{n-y}$$

于是我们得到了后验分布

$$p(\theta\lvert y)\propto \theta^y (1-\theta)^{n-y}$$

即 $\theta\lvert y\sim\mathrm{Beta}(y+1,n-y+1)$

如果读者对于$\mathrm{Beta}$分布并不熟悉,请参看Beta distribution

下图展示了当数据规模变大而为$1$的比例都为0.6时,后验分布的变化情况。

绘图代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
par(mfrow=c(2,2))
plot_n <- function(n) {
y = n * 0.6
f <- function(x) { dbeta(x, y+1, n-y+1)}
curve(f,0,1)
# the code below is just for plot
maximum = dbeta(y/n,y+1, n-y+1)
text(x=0.8, y=0.8 * maximum, labels=paste0('n=',3))
text(x=0.8, y=0.6 * maximum, labels=paste0('y=',y))
}
plot_n(5)
plot_n(20)
plot_n(100)
plot_n(1000)

可以看到随着数据的增大,我们的后验分布更加集中在数据里给出的信息$0.6$。

此时,我们使用后验分布的均值作为点估计,有

$$\mathbb{E}(\theta\lvert y)=\frac{\alpha}{\alpha+\beta}=\frac{y+1}{n+2}$$

事实上,对于后验均值有这样的看待方式

$$\frac{y+1}{n+2}=\frac{n}{n+2}\frac{y}{n}+\frac{2}{n+2}\frac{1}{2}$$

其中$y/n$,$1/2$分别是数据和先验给出的均值,而后验均值则是二者使用样本量大小所做的一个加权平均。

为什么先验分布的样本量会是$2$呢?因为均匀分布其实是$\mathrm{Beta}(1,1)$,相当于自带两个数据且一正一负,所以后验分布会成为那样加权的形式。

考虑我们之前对于传统统计方法的批评——过于依赖数据的结果。而先验的引入正是为了将我们从数据“拉开”,而同时考虑先验知识的影响。甚至在共轭先验下,先验通常可以解释为额外的数据(后面会介绍)。在开头的例子中,2次实验得到2次1后,我们的后验均值为$3/4$,而20000次实验得到20000次1的后验均值为$20001/20002$,这样的结果看起来就更加自然了:随着实验的增加,当我发现出现的结果都是1时,我心中对于参数的认识也就越来越趋近于$1$。值得注意的是,在这样的框架下后验均值永远不会是$1$,这也是可以理解的——我就目前的数据来说倾向于参数为$1$,但由于仍有未观测到的数据,我仍然给其留了一份余地。

后验分布的性质

使用概率中的全期望公式Law of total expectation和条件方差公式,我们有如下关系

$$\mathbb{E}(\theta)=\mathbb{E}(\mathbb{E}(\theta\lvert y))$$
$$\mathrm{Var}(\theta)=\mathbb{E}(\mathrm{Var}(\theta\lvert y))+\mathrm{Var}(\mathbb{E}(\theta\lvert y))$$

这里不想从数学的角度证明与介绍这些式子,让我们在贝叶斯框架审视这两条公式。

  • 第一条说明:先验分布的均值$\mathbb{E}(\theta)$是所有后验分布的均值$\mathbb{E}(\theta\lvert y)$在数据分布$p(y)$上的加权平均。

  • 第二条说明:后验分布的方差$\mathrm{Var}(\theta\lvert y)$在平均意义上(数据分布$p(y)$上的加权平均)不会超过先验方差$\mathrm{Var}(\theta)$。

    • 其直观理解是,输入了信息不会使得我们的预测的方差更大
    • 值得注意的是,其只保证了平均意义上后验方差不会更大。对于某一具体的观测$y$,其后验方差是可以超过先验方差的。
  • 第二条还说明:后验均值的方差越大,$\theta$的方差减少(对比先后验)就更多。

先验的选择

在本文的引子里,我们介绍了二项分布的例子。现在让我们选用另一个先验分布$\theta\sim\mathrm{Beta}(\alpha,\beta)$,可以推导出其后验分布$\theta\lvert y\sim\mathrm{Beta}(\alpha,\beta)$。可以看到先验的变化会影响后验的分布以及与之相关的推断,因此选取一个合适的先验就十分重要了。其中,共轭先验和无信息先验是广泛被应用的两种方法。

共轭先验

定义

共轭先验的定义如下:如果$\mathcal{F}$是一个分布族,$\mathcal{P}$是这个分布的参数的分布族。如果对于所有$\mathcal{F}$中的概率密度函数$p(\cdot\lvert\theta)$和$\mathcal{P}$中的先验$p(\theta)$,我们都有后验分布$p(\theta\lvert y)\in\mathcal{P}$,我们就称$\mathcal{P}$对$\mathcal{F}$是共轭的。简而言之,先验和后验属于同一个函数族。

优劣

共轭先验的好处如下

  • 先验和后验有相同的形式:因为均属于同一个分布族。这通常意味着我们能给出一个闭式解
  • 先验和后验有相同的形式,意味着我们可以将后验作为下一时刻的先验:因此我们可以不断根据新来的数据然后调整后验分布的参数,也就是所谓的在线贝叶斯学习(online bayesian learning)
  • 通常,我们可以将先验解释为额外的数据,这使得我们的模型有良好的解释性

既然共轭先验这么好,那我们为什么还需要非共轭的先验呢?

  • 实际应用中,复杂模型可能无法求出共轭先验
  • 使用非共轭先验并没有带来任何的问题
  • 许多非共轭先验其实是使用共轭先验的混合构造出来的,因为单纯的共轭先验可能并不合理

那我们如何求解共轭先验呢?幸运的是,对于指数族分布,其共轭分布一定存在。(如果读者对指数族分布、自然参数、充分统计量等名词并不熟悉,请参考我关于统计推断的笔记(Not finished yet so use google instead now).

求法

对于指数族分布

$$p(y_i\lvert\theta)=f(y_i)g(\theta)\exp\{\phi(\theta)^\intercal u(y_i)\}$$

我们可以得到概率密度函数

$$p(y\lvert\theta)=\left(\prod_{i=1}^{n}f(y_i)\right)g(\theta)^n\exp\{\phi(\theta)^\intercal \sum_{i=1}^nu(y_i)\}\propto g(\theta)^n\exp\{\phi(\theta)^\intercal t(y)\}$$

而我们的先验分布取为$p(\theta)\propto g(\theta)^\eta\exp\{\phi(\theta)^\intercal \nu$,则后验分布为

$$p(\theta\lvert y)\propto g(\theta)^{\eta+n}\exp\{\phi(\theta)^\intercal (\nu+t(y))\}$$

当然,通常我们只需观察概率密度函数的形式,就能写出共轭先验的形式,具体可以参见下面的例子。

例子

已知方差的正态分布

先考虑只有一个样本的简单情形。由于概率密度函数有形式$p(y\lvert\theta)\propto\exp\{-(y-\theta)^2/(2\sigma^2)\}$

我们取先验函数$p(\theta)\propto \exp\{A\theta^2+B\theta+C\}$,或等价的,$\exp\{-(\theta-\mu_0)^2/(2\tau_0^2)\}$

可以算出后验分布

$$p(\theta\lvert y)\propto\exp\left\{-\frac{1}{2}\left(\frac{(y-\theta)^2}{\sigma^2}+\frac{(\theta-\mu_0)^2}{\tau_0^2}\right)\right\}$$

整理得$\theta\lvert y\sim\mathcal{N}(\mu_1,\tau_1^2)$,其中

$$\mu_1=\frac{\frac{\mu_0}{\tau_0^2}+\frac{y}{\sigma^2}}{\frac{1}{\tau_0^2}+\frac{1}{\sigma^2}}\quad\mathrm{and}\quad\frac{1}{\tau_1^2}=\frac{1}{\tau_0^2}+\frac{1}{\sigma^2}$$

可以看到熟悉的,先验和后验的加权平均的形式。

让我们看一下后验预测分布(posterior predictive distribution)

$$p(\tilde{y}\lvert y)=\int p(\tilde{y}\lvert \theta)p(\theta\lvert y)\,\mathrm{d}\theta\propto\int \exp\left\{-\frac{(\tilde{y}-\theta)^2}{2\sigma^2}\right\}\exp\left\{-\frac{(\theta-\mu_1)^2}{2\tau_1^2}\right\}\,\mathrm{d}\theta$$

形式比较复杂,但求期望和方差的时候无需使用该密度函数,因为

$$\mathbb{E}(\tilde{y}\lvert y)=\mathbb{E}(\mathbb{E}(\tilde{y}\lvert \theta,y)\lvert y)=\mathbb{E}(\theta\lvert y)=\mu_1$$

$$\mathrm{Var}(\tilde{y}\lvert y)=\mathbb{E}(\mathrm{Var}(\tilde{y}\lvert y)\lvert y)+\mathrm{Var}(\mathbb{E}(\theta\lvert y)\lvert y)=\mathbb{E}(\sigma^2\lvert y)+\mathrm{Var}(\theta\lvert y)=\sigma^2+\tau_1^2$$

多样本情形也完全没有变复杂,有$\theta\lvert y_1,\ldots,y_n\sim\mathcal{N}(\mu_n,\tau_n^2)$,其中

$$\mu_n=\frac{\frac{1}{\tau_0^2}\mu_0+\frac{n}{\sigma^2}\bar{y}}{\frac{1}{\tau_0^2}+\frac{n}{\sigma^2}}\quad\mathrm{and}\quad\frac{1}{\tau_n^2}=\frac{1}{\tau_0^2}+\frac{n}{\sigma^2}$$

具体的计算细节留给读者自行解决。当然,如果发现可以将$y_1,\ldots,y_n$使用充分统计量$\bar{y}$替代的话,其结论是显然的。

已知均值的正态分布

直接上多样本的情形了。首先写出概率密度函数

$$p(y\lvert\sigma^2)\propto \sigma^{-n}\exp\left\{-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\theta)^2\right\}\overset{\Delta}{=}(\sigma^2)^{-n/2}\exp\left\{\frac{-n\nu}{2\sigma^2}\right\}$$

根据这个形式,我们取共轭先验逆伽马分布Inverse-Gamma distribution,取$\alpha=\nu_0/2,\beta=\nu_0\sigma_0^2/2$

一个等价的写法是带尺度的逆卡方分布Inverse-chi-squared distribution,其关系式子为$\mathrm{Inv-}\chi^2(\nu,s^2)\equiv\mathrm{Inv-Gamma}(\nu/2,\nu s^2/2)$,此时先验的形式更简洁——$\sigma^2\sim \mathrm{Inv-}\chi^2(\nu_0,\sigma_0^2)$。本书和本笔记都采取了这个记号。对于初学者可能对此有点难以接受,但习惯了后这样的记号还是十分方便的。

此时先验分布的形式为

$$p(\sigma^2)\propto(\sigma^2)^{-(\nu_0/2+1)}\exp\left\{-\frac{\nu_0\sigma_0^2}{2\sigma^2}\right\}$$

$$\begin{align}
p(\sigma^2\lvert y)
& \propto p(\sigma^2)p(y\lvert\sigma^2)\\
& \propto (\sigma^2)^{-(\nu_0/2+1)}\exp\left\{-\frac{\nu_0\sigma_0^2}{2\sigma^2}\right\}(\sigma^2)^{-n/2}\exp\left\{\frac{-n\nu}{2\sigma^2}\right\}\\
& \propto (\sigma^2)^{-((n+\nu_0)/2+1)}\exp\left\{-\frac{\nu_0\sigma_0^2+n\nu}{2\sigma^2}\right\}
\end{align}$$

即$\sigma^2\lvert y\sim \mathrm{Inv-}\chi^2((n+\nu_0),(\nu_0\sigma_0^2+n\nu)/(\nu_0+n))$

或等价的,$\mathrm{Inv-Gamma}((n+\nu_0)/2,(\nu_0\sigma_0^2+n\nu)/2)$

泊松分布

概率密度函数

$$p(y\lvert\sigma^2)=\prod_{i=1}^{n}\frac{\theta^{y_i}}{y_i!}\exp\{-\theta\}\propto \theta^{t(y)}\exp\{n\theta\},\,t(y)=\sum_{i=1}^n y_i=n\bar{y}$$

于是取共轭先验为伽马分布Gamma distribution

$$p(\theta)\propto \theta^{\alpha-1}\exp\{\beta\theta\}$$

由此可以求得后验分布

$$p(\theta)\propto \theta^{\alpha+n\bar{y}-1}\exp\{\beta\theta+n\theta\},\,\theta\lvert y\sim\mathrm{Gamma}(\alpha+n\bar{y},\beta+n)$$

共轭先验的应用

当选定了共轭先验后,剩下的就是确定先验分布参数的事情了。通常地,如果你之前有相关的数据,一个常用的方法就是使用矩估计,将样本矩使用超参数表示出来,最后通过解方程(组)确定出超参数的信息。如果没有这样的数据信息,也可以通过询问相关的专家一系列问题,比如峰值大概会在哪里取得,方差大致有多大等信息来约束出一个合理的超参数。当然,先验的准确性通常对结果的影响并没有特别大,尤其当数据量较大的时候,因此实际中也无需太过纠结这一点。

无信息先验

如上所言,但当我们既没有先前的数据,也没有专家提供信息的时候,我们通常就会选取无信息先验:我们不对结果做任何主观的偏好(当然,后面我们会看着这个说法还有待商榷),而让数据自己说话(let the data speak for themselves)。是的,我想许多读者也会猜到——均匀分布算得上是一个特别“无信息”的先验了。但无信息先验远不止这么简单。

合适的先验

均匀分布有一个最大的问题,就是我们无法在一个无限区间上定义均匀分布。因此,对于参数空间是无限区间的情况,使用均匀分布作为无信息先验看似是不可行的。然而,很多时候这种情况是可以被处理的,为此我们先引入“合适”这个概念。

简单说,一个合适(proper)的先验意味着其是一个合理的概率分布,简而言之就是满足非负性且其全空间积分为1(或有限,因为有限值总可以使用一个归一化因子将其化为1)。但是对于一个Improper的先验,其后验是有可能proper的,在这样的情况下我们也是可以使用这样的先验的。当然,这样的关系不一定成立。

回忆上面谈到的已知方差的正态分布,我们有$\theta\lvert y\sim\mathcal{N}(\mu_1,\tau_1^2)$,其中

$$\mu_1=\frac{\frac{\mu_0}{\tau_0^2}+\frac{y}{\sigma^2}}{\frac{1}{\tau_0^2}+\frac{1}{\sigma^2}}\quad\mathrm{and}\quad\frac{1}{\tau_1^2}=\frac{1}{\tau_0^2}+\frac{1}{\sigma^2}$$

我们可以将均匀分布看为正态分布的极限——当正态分布的方差足够大的时候,整个曲线好似被拉平并压扁一般,也就形成了“均匀分布”,因此在上式中我们将$\tau_0^2\rightarrow\infty$,我们立得$\theta\lvert y\sim\mathcal{N}(\bar{y},\sigma^2/n)$

读者应当对这一结果并不意外:这与传统统计学里面得到的结果一致,因为此时只有数据在说话,先验的影响被我们压下去了。

当然,事实上我们并不需要走这样一个复杂的推导过程再取极限,我们先取先验分布为$f(\theta)=1$,显然这是一个Improper的先验因为其全空间积分不存在。然而,将其带进去算后验分布得到的结果是proper的——正是上面的结果。因此我们可以选择这样的一个improper的平坦先验作为我们的先验分布。

均匀先验的问题

看似均匀先验基本解决了无信息先验的问题,但我们需要注意到如下事实

设$\theta$是一个随机变量并有概率密度函数$p(\theta)$,记$\phi=h(\theta)$是一个一一映射,此时我们去计算$\phi$的密度函数,有

$$f(\phi)=p(\theta)\Bigg\lvert\frac{\mathrm{d}\theta}{\mathrm{d}\phi}\Bigg\rvert=p(\theta)\lvert h^\prime(\theta)\rvert^{-1}$$

这说明,如果$\theta$是均匀分布,那么$\phi$通常就不再是均匀分布了,一个具体的例子是,如果$\sigma$服从均匀先验,那么$p(\sigma^2)=(2\sigma)^{-1}$就不再是均匀分布,而且可以看出这样的先验更加青睐小方差。同样的,如果我们让$\sigma^2$服从均匀先验,那么$p(\sigma)=2\sigma$,青睐大方差。

因此,事实上不存在一个一致的均匀先验,为此人们也提出过一些对应的解决方法。

Jeffrey 先验

沿用上面的记号,并引入一个新的概念:Fisher信息量(如果不了解,请参考我关于统计推断的笔记(Not finished yet so use google instead now))

$$J(\theta)=\mathbb{E}\left(\left(\frac{\mathrm{d}\log p(y\lvert\theta)}{\mathrm{d}\theta}\right)\big\rvert\theta\right)=-\mathbb{E}\left(\frac{\mathrm{d}^2\,\log p(y\lvert\theta)}{\mathrm{d}\theta^2}\big\rvert\theta\right)$$

我们同时求取$J(\phi)$

$$J(\phi)=-\mathbb{E}\left(\frac{\mathrm{d}^2\,\log p(y\lvert\phi)}{\mathrm{d}\phi^2}\Bigg\rvert\phi\right)=-\mathbb{E}\left(\frac{\mathrm{d}^2\,\log p(y\lvert|theta=h^{-1}(\phi))}{\mathrm{d}\theta^2}\Bigg\lvert\frac{\mathrm{d}\theta}{\mathrm{d}\phi}\Bigg\rvert^2\Bigg\rvert\theta\right)=J(\theta)\Bigg\lvert\frac{\mathrm{d}\theta}{\mathrm{d}\phi}\Bigg\rvert^2$$

即$$J(\phi)^{1/2}=J(\theta)^{1/2}\Bigg\lvert\frac{\mathrm{d}\theta}{\mathrm{d}\phi}\Bigg\rvert$$

于是我们可以选择先验$p(\theta)\propto J(\theta)^{1/2}$,此时$p(\phi)=p(\theta)\lvert h^\prime(\theta)\rvert^{-1}\propto J(\phi)$。可以看到此时先验的形式不会随着参数形式的变换而变换。

以上文提到的方差为例,使用Jeffrey先验可以得到

$$p(\sigma)\propto \frac{1}{\sigma},\qquad p(\sigma^2)\propto \frac{1}{\sigma^2},\qquad p(\log \sigma^2)\propto p(\log \sigma)\propto 1$$

另一个例子就是最早提到的二项分布的例子

  • 使用无信息先验,则$p(\theta)=1$,即$\theta\sim\mathrm{Beta}(1,1)$
  • 使用Jeffrey先验,则$p(\theta)\propto\theta^{-1/2}(1-\theta)^{-1/2}$,即$\theta\sim\mathrm{Beta}(1/2,1/2)$
  • 使用自然参数的无信息先验,则$p(\mathrm{logit}(\theta))=1$,即$\theta\sim\mathrm{Beta}(0,0)$
    • 如果不了解,请参考我关于统计推断的笔记(Not finished yet so use google instead now),但这不重要

位置族参数和尺度族参数

以位置族参数为例,如果$X$密度函数可以写成关于$f(x-\theta)$,则称$\theta$为位置族参数($\theta\in\mathbb{R}$)。考察下列变换$Y=X+c,\quad \eta=\theta+c$,此时$Y$的密度函数为$p(y-\eta)$。

一方面,我们希望$\theta$和$\eta$的先验具有相同的形式,即$\pi(\eta)=\pi^{\ast}(\eta)$

另一方面由于有关系$\eta=\theta+c$,我们有$\pi^{\ast}(\eta)=\lvert\frac{\mathrm{d}\theta}{\mathrm{d}\eta}\rvert\pi(\theta)=\pi(\theta)=\pi(\eta-c)$

因此$\pi(\eta)=\pi(\eta-c)$。而$\eta,c$都可以任意选择,也就意味着$\pi(\eta)=\mathrm{constant}$,也就是均匀分布。

尺度族参数的推导完全类似,其定义为如果$X$密度函数可以写成关于$\frac{1}{\sigma}f(\frac{x}{\sigma})$,则称$\sigma$为尺度族参数($\sigma>0$)

同样的,我们有约束

$\pi(\sigma)=\pi^{\ast}(\sigma)$,且$\pi^{ast}(\sigma)=\lvert\frac{\mathrm{d}\theta}{\mathrm{d}\sigma}\rvert\pi(\theta)=\frac{1}{c}\pi(\frac{\sigma}{c})$

于是$\pi(\sigma)=\frac{1}{c}\pi(\frac{\sigma}{c})$。不妨取$c=1,\pi(1)=1$,则$\pi(\sigma)=1/\sigma$,也就是Jeffrey先验。

通常,对于这两族参数我们优先考虑这两个先验。其他的参数我们可能根据需要确定是使用无信息还是Jeffery,甚至有些时候需要根据后验分布(为了使后验分布proper)而选取合理的先验。