3.3 参数估计：贝叶斯估计

3.3 贝叶斯估计

矩估计和极大似然估计方法的优点是比较客观客观，基本由随机采样数据决定。缺点是需要在大样本情况下估计才比较准确。不能把人类知识用于估计。例如，某公司研发新产品，需要估计合格率，这是典型的伯努利分布。按照矩估计和极大似然估计方法，需要试生产大量产品后才能获得比较好的估计，这在实践中十分昂贵和耗时。该公司的研发人员根据同类产品的历史经验和理论分析或仿真，可以对新产品的合格率有个比较可靠的估计，即使新产品还没有生产。这些知识是十分宝贵的，如果不能充分利用则实在可惜。如果很好的利用了这些知识，则可能只需少量样本即可获得很好的估计，贝叶斯估计就是这样一种方法。

贝叶斯估计有个核心概念：先验。先验是先于试验即在试验之前，注意不是实验。在统计估计中，试验指通过对总体分布进行不断抽样，根据抽样来获得对总体分布的知识。例如要估计黑箱中黑球白球的比例，通过不断取球，根据取出的球中黑球白球比例来估计真实比例。总之抽样获得的知识就是经验，所以在统计估计中先验指在抽样之前。先验知识指在抽样之前就存在的关于总体分布的知识：总体分布和分布参数，这些知识如何获得的？有可能是前人研究类似分布获得的知识，也有可能是通过理论分析获得的知识等等。还是以估计新产品合格率为例，先验知识就是研发人员根据同类产品和经验，或者对新产品进行理论研究获得的关于合格率的知识。注意人们关于这个合格率，其实并不是一个固定的值，而是一个分布。例如如果根据先验知识，知道新产品合格率会比较高，则会认为合格率为0.8的概率会远大于0.1的概率，如合格率在 $[0.8,1]$ 的概率为0.7，在 $[0,0.1]$ 的概率为0.1。在抽样之前，关于合格率分布的密度函数就是先验分布，这是贝叶斯估计的核心概念。矩估计和极大似然估计都是假设合格率是一个固定值而不是分布，只是我们不知道这个真实值而已。是否利用先验知识和被估计值是否是概率分布，这是贝叶斯估计和矩估计和极大似然估计最大区别，为此贝叶斯估计被称为贝叶斯学派，矩估计和极大似然估计等称为频率学派。两个学派各有优缺点，历史上这两种方法争论不休。

知道了先验分布就很容易理解贝叶斯估计了。根据条件概率公式

$$p(A|B)=p(A,B)/p(B)$$

$p(A|B)$ 表示事件B发生条件下事件A发生的概率，$p(A,B)$ 表示事件A，B同时发生的概率。同理可得

$$p(B|A)=p(A,B)/p(A)$$

故 $p(A|B)p(B)=p(B|A)p(A)$ 即

$$p(B|A)=p(B)[p(A|B)/p(A)]$$

这就是著名的贝叶斯公式，贝叶斯估计就是依据该公式。

怎么理解贝叶斯公式呢？$p(B)$ 是事件B发生概率，$p(B|A)$ 是事件A发生条件下事件B发生概率，他们之间的关系由 $p(A|B)/p(A)$ 决定，这个表达式体现了事件A对事件B的影响。如果 $p(A|B)/p(A)=1$ 则 $p(B|A)=p(B)$ ，这说明事件A对事件B没有影响，他们是独立的。根据 $p(A|B)/p(A)=1$ 得 $p(A,B)=p(A)p(B)$ ,A,B事件确实是独立的。如果 $p(A|B)/p(A)>1$ 则 $p(B|A)>p(B)$ ，这说明事件A对事件B的发生有促进作用。由 $p(A|B)/p(A)>1$ 得 $p(A|B)>p(A)$ ，这说明事件B对事件A的发生有促进作用，所以事件A，B是相互促进的，和生活常识相符。同理如果 $p(A|B)/p(A)<1$ 则 $p(B|A)<p(B)$ ，这说明事件A对事件B的发生有抵消作用。由 $p(A|B)/p(A)<1$ 得 $p(A|B)<p(A)$ ，这说明事件B对事件A的发生有抵消作用，所以事件A，B是相互抵消的。

贝叶斯公式中事件A，B可以是任意事件，对于贝叶斯估计来说，事件A就是对总体随机抽样获得 $n$ 个样本 $A:{\mathbf{x}}_{data}=(x_1,\cdots ,x_n)$ ，是经验。事件B就是要估计的参数向量 $B=\theta$ 。$p(B)=p(\theta )$ 是参数向量的概率分布，是未试验之前给出的，就是先验分布，在贝叶斯估计中记为 $\pi (\theta )$ 。$p(B|A)=p(\theta |{\mathbf{x}}_{data})$ 是采样数据后参数向量的概率分布，称为后验分布，在贝叶斯估计中记为 $\pi (\theta |{\mathbf{x}}_{data})$ 。后验分布就是有了数据后对参数向量的先验分布进行更新后得到的，包含了先验知识和采样信息。$p(A|B)=p({\mathbf{x}}_{data}|\theta )$ 是假设参数向量 $\theta$ 已知情况下采样到数据 ${\mathbf{x}}_{data}$ 的概率，就是极大似然估计方法中的似然函数！这点十分重要，建立了和极大似然估计的关系。$p(A)=p({\mathbf{x}}_{data})$ 就是采样到数据 ${\mathbf{x}}_{data}$ 的边缘概率。所以贝叶斯估计的公式为

$$\pi (\theta |{\mathbf{x}}_{data})=\frac{\pi (\theta )p({\mathbf{x}}_{data}|\theta )}{p({\mathbf{x}}_{data})}$$

对连续型随机变量有

$$p({\mathbf{x}}_{data})={\int }_{\theta }\pi (\theta )p({\mathbf{x}}_{data}|\theta )d\theta$$

贝叶斯估计公式看起来很简单，但实际求解非常复杂。一般来说，计算 $p({\mathbf{x}}_{data})$ 是不太可能的，因为需要计算积分，需要精心选择合适的先验分布才有可能计算积分。

对于伯努利分布来说，需要估计随机变量取 1 的概率 $p$ ，如果 $p=\theta$ 的先验分布假设为Beta分布，Beta分布的概率密度函数为：

$$beta(\theta ;(\alpha ,\beta ))=\frac{1}{B(\alpha ,\beta )}{\theta }^{\alpha -1}(1-\theta {)}^{\beta -1};0<\theta <1$$

其中 $\alpha ,\beta$ 是分布参数，需要我们利用先验知识获得 $\theta$ 的分布，然后选择最优的 $\alpha ,\beta$ 来拟合该分布。

Beta分布的两个重要性质，分布众数（函数极值点）为 $\frac{\alpha -1}{\alpha +\beta -2}$ ，分布均值为 $\frac{\alpha }{\alpha +\beta }$ ；$beta(\theta ;(\alpha =1,\beta =1))=1$ 是均匀分布。

假设随机抽样 $n$ 个样本，$m$ 次随机变量取 1，则似然函数为

$$p({\mathbf{x}}_{data}|\theta )={\theta }^m(1-\theta {)}^{n-m}$$

抽样数据的边缘密度为

$$\begin{gathered} p({\mathbf{x}}_{data})={\int }_{\theta }\pi (\theta )p({\mathbf{x}}_{data}|\theta )d\theta \\ ={\int }_{\theta }\frac{1}{B(\alpha ,\beta )}{\theta }^{\alpha -1}(1-\theta {)}^{\beta -1}{\theta }^m(1-\theta {)}^{n-m}d\theta \\ =\frac{B(\alpha +m,\beta +n-m)}{B(\alpha ,\beta )} \end{gathered}$$

注意边缘分布与参数 $\theta$ 无关！

所以参数的后验分布为

$$\begin{gathered} \pi (\theta |{\mathbf{x}}_{data})=\frac{\pi (\theta )p({\mathbf{x}}_{data}|\theta )}{p({\mathbf{x}}_{data})} \\ =\frac{\frac{1}{B(\alpha ,\beta )}{\theta }^{\alpha -1}(1-\theta {)}^{\beta -1}{\theta }^m(1-\theta {)}^{n-m}}{\frac{B(\alpha +m,\beta +n-m)}{B(\alpha ,\beta )}} \\ =beta(\theta ;(\alpha +m,\beta +n-m)) \end{gathered}$$

后验分布也是Beta分布，和先验分布属于同类分布。当后验分布和先验分布属于同族分布时，此时称先验分布为共轭先验分布。采用共轭先验分布可以简化计算。

有了后验分布就可以进行各种估计了，可以进行点估计和区间估计。

点估计有后验分布的众数（后验密度最大点）、后验分布的中位数、后验分布的均值。其中后验均值的计算公式为 ${\int }_{\theta }\theta \pi (\theta |{\mathbf{x}}_{data})d\theta$ 。那么这三种估计哪一种最好呢？各有优缺点，实践中都可以采用，后验分布的均值在平均意义下是最优的。如果后验分布是中心对称且中心点是极大值，则这三种点估计是同一值。

区间估计可以计算参数在任意区间的概率 $p={\int }_{\theta \in \Omega }\theta \pi (\theta |{\mathbf{x}}_{data})d\theta$ ，很可惜由于后验概率很复杂，该积分式一般没有解析值，只能数值计算或查表。历史上贝叶斯本人可能就是因为没有计算出Beta分布区间概率的准确值，在生前迟迟没有公开发表贝叶斯估计论文，直到死后两年才发表。或者计算最高后验概率密度区间(highest posterior density HPD)，即希望在置信水平 $1-\alpha$ 下，置信区间的长度越短越好。当参数的后验密度函数 $\pi (\theta |{\mathbf{x}}_{data})$ 曲线是单峰对称的，则峰值两侧各取 $\alpha /2$时，置信区间的长度为最短，此置信区间为 HPD 置信区间。当后验密度曲线为单峰非对称时， 则 HPD 置信区间求起来相对复杂。最高后验概率密度区间的意义可以这样理解，被估计参数的真值有 $1-\alpha$ 概率落在该区间，虽然贝叶斯估计认为被估计值不存在真值，只存在一个分布，但这种理解非常符合人的直觉。

贝叶斯估计从一定意义上说也遵循大数定理和中心极限定理。即随着抽样样本数量的增加，后验分布函数曲线越来越集中，只形成一个尖峰，后验分布方差越来越小，当样本无穷多时，方差趋近 0 。而且，随着抽样样本数量的增加，先验分布对后验分布的影响越来越小，当样本无穷多时，先验分布影响趋于 0。这个十分符合直观，因为样本无穷多时，根据大数定理和中心极限定理，样本信息足够估计参数，完全不需要先验知识，先验知识只在样本比较少的时候才有用，样本越少，先验知识越重要，样本越多先验知识越不重要，毕竟先验知识主观成分比较大，抽样样本数据更有说服力。样本无穷多时分布极值点必然和极大似然估计结果一致。

还是以伯努利分布为例，后验分布为

$$\pi (\theta |{\mathbf{x}}_{data})=beta(\theta ;(\alpha +m,\beta +n-m))$$

后验分布均值为 $E(\theta )=\frac{\alpha +m}{\alpha +\beta +n}$，分布众数为 $\frac{\alpha +m-1}{\alpha +\beta +n-2}$，后验方差为 $Var(\theta )=\frac{(\alpha +m)(\beta +n-m)}{(\alpha +\beta +n{)}^2(\alpha +\beta +n+1)}$ 。

后验分布均值可变换为

$$\begin{gathered} E(\theta )=a\frac{m}{n}+(1-a)\frac{\alpha }{\alpha +\beta } \\ a=\frac{n}{\alpha +\beta +n} \end{gathered}$$

可见后验分布均值是样本均值 $\frac{m}{n}$ 和先验均值 $\frac{\alpha }{\alpha +\beta }$ 的加权平均，权重为 $a$ ，可见后验均值综合了样本信息和先验信息。显然当 $n$ 越来越大时权重 $a$ 越来越接近 1，后验均值越来越接近样本均值，样本均值就是极大似然估计值，与先验分布无关。当样本数量少时，权重 $a$ 小，先验知识很重要。从另一个角度看先验分布 $beta(\theta ;(\alpha ,\beta ))$ ，可以认为是进行了虚拟试验，总共试验了 $\alpha +\beta$ 次，其中 $\alpha$ 次取值为1 即成功次数。

后验方差可变换为

$$Var(\theta )=\frac{E(\theta )(1-E(\theta ))}{\alpha +\beta +n+1}$$

显然当 $n$ 越来越大时方差越来越小，无穷多时，方差趋近 0，这表示估计参数值几乎收敛到一个点，失去了随机性，估计值是百分之百准确的。 $n$ 小时，方差较大，估计准确度不高。

当样本分布为高斯分布时 $N(\mu ,{\sigma }_0^2)$ ，其中 $\mu$ 未知，${\sigma }_0^2$ 已知。若 $\mu$ 的先验分布为高斯分布 $N({\mu }_b,{\sigma }_b^2)$ 。根据贝叶斯估计公式，可以计算得到 $\mu$ 的后验分布也是高斯分布，$N({\mu }_p,{\sigma }_p^2)$ 。

$$\begin{gathered} {\mu }_p=\frac{\bar{x}{h}_0+{\mu }_b{h}_1}{h_0+h_1} \\ ({\sigma }_p^2{)}^{-1}=h_0+h_1 \\ {h}_0=({\sigma }_0^2/n{)}^{-1} \\ {h}_1=({\sigma }_b^2{)}^{-1} \end{gathered}$$

可见，在方差已知情况下高斯分布的共轭先验分布为高斯分布。后验均值是样本均值 $\bar{x}$ 和先验均值的加权平均，后验方差是样本平均方差和先验方差的调和平均。显然当 $n$ 越来越大时 $h_0$ 越来越大，后验均值越来越接近样本均值，样本均值就是极大似然估计值，与先验分布无关。当样本数量少时，先验知识很重要。当 $n$ 越来越大时方差越来越小，无穷多时，方差趋近 0，这表示估计参数值几乎收敛到一个点，失去了随机性，估计值是百分之百准确的。 $n$ 小时，方差较大，估计准确度不高。

关于点估计中后验分布众数，也称为极大后验估计MAP。根据后验公式

$$\pi (\theta |{\mathbf{x}}_{data})=\frac{\pi (\theta )p({\mathbf{x}}_{data}|\theta )}{p({\mathbf{x}}_{data})}$$

由于数据边缘分布 $p({\mathbf{x}}_{data})$ 是对参数积分所得，与参数 $\theta$ 无关。故后验分布最大值等价于 $\pi (\theta )p({\mathbf{x}}_{data}|\theta )$ 最大值，是似然函数和先验分布的乘积。两边取对数得

$$logp({\mathbf{x}}_{data}|\theta )+log\pi (\theta )$$

故极大后验估计与极大似然估计十分相似，只是多了先验分布对数。当先验分布是常数分布时，极大后验估计与极大似然估计结果一致。从广义上看，$log\pi (\theta )$ 是正则项。似然函数 $logp({\mathbf{x}}_{data}|\theta )$ 是样本信息，$log\pi (\theta )$ 是先验知识。

现在分析贝叶斯学派和频率学派的优缺点。进行参数估计，贝叶斯学派利用了三种信息：先验知识获得先验分布；采样的样本信息；和一个不太明显的信息：样本总体分布假设，比如假设样本分布满足高斯分布或拉普拉斯分布，而不是其他分布，这也是需要先验知识的。频率学派利用了其中两种信息：采样的样本信息；样本总体分布假设。

贝叶斯学派最大争议是先验分布，由于选择先验分布具有主观性，表现出一定的随意性，导致分析失去部分客观性。不同的先验分布对后验分布有影响，甚至有时候，两个『接近的』先验分布会导致很不同的后验分布。所以确定合理的先验分布是贝叶斯学派的核心。如果先验知识能定出合理的先验分布，此时贝叶斯估计比频率学派能获得更好结果，特别是样本数量少时。当先验很弱或样本数量很大，不需要勉强采用贝叶斯估计。后验分布难以计算也是贝叶斯实际应用中的拦路虎。贝叶斯学派的优点是利用后验分布进行参数估计十分符合人们的直觉。

频率学派的缺点是完全忽略先验知识，当有很多先验知识时这会显得十分不明智，特别是样本很少时。频率学派进行参数估计时难以理解，不太符合人们直觉。优点是客观性比较高，但选择合适的样本总体分布假设，实际上也有一定的主观性，所以不存在完全客观的统计推断方法。

计算后验分布的难题，现在采用马尔可夫链蒙特卡洛模拟 MCMC 可部分解决，具体内容读者可查阅相关教材。

关于先验分布的选择，我们再强调一下。当我们不存在任何先验知识时，但还是希望用贝叶斯估计。那怎么确定先验分布呢？此时先验分布称为无信息先验分布。

还是以伯努利分布为例，我们对参数 $p$ 一无所知，显然最合理的先验分布为均匀分布，即 $(0,1)$ 内均匀分布，$beta(\theta |(\alpha =1,\beta =1))$ 就是均匀分布。利用该分布可以获得合理的后验分布，比极大似然估计更合理，相当于在正式试验前进行了两次虚拟试验，其中成功了一次。

高斯分布的均值的先验分布是什么呢？由于均值可以取任意实数，所以一个合理的先验分布是取整个实数域上的均匀分布，显然此时先验分布在定义域内积分为无穷大，不满足先密度函数在定义域内积分为 1 的要求。严格来说不能作为先验分布。但实在找不出其他更合理的先验分布，怎么办？幸好，采用整个实数域上的均匀分布计算的后验分布居然是一个真正概率分布，我们真正关心的是后验分布，所以我们此时定义这种先验分布为广义先验分布。

更有趣的是，高斯分布方差 ${\sigma }^2$ 的平方根标准差 $\sigma$ 先验分布是什么呢？也是 $(0,\infty )$ 的均匀分布吗？其实不是，$\pi (\sigma )={\sigma }^{-1}$ 更合理，神奇吧！关于各种设置先验分布的知识读者可查阅相关教材。