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 | par(mfrow=c(2,2)) |
可以看到随着数据的增大，我们的后验分布更加集中在数据里给出的信息\(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）而选取合理的先验。