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

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

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

模型假设

我们从二项分布中抽取了一系列样本\(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)而选取合理的先验。