6.5 贝叶斯估计

在统计学中有两个大的学派:频率学派(也称经典学派)和贝叶斯学派.本书主要介绍频率学派的理论和方法,此一小节将对贝叶斯学派做些介绍.

6.5.1 统计推断的基础

我们在前面已经讲过,统计推断是根据样本信息对总体分布或总体的特征数进行推断,事实上,这是经典学派对统计推断的规定,这里的统计推断使用到两种信息:总体信息和样本信息;而贝叶斯学派认为,除了上述两种信息以外,统计推断还应该使用第三种信息:先验信息.下面我们先把三种信息加以说明.

1. 总体信息

总体信息即总体分布或总体所属分布族提供的信息.譬如,若已知“总体是正态分布”,则我们就知道很多信息.譬如:总体的一切阶矩都存在,总体密度函数关于均值对称,总体的所有性质由其一、二阶矩决定,有许多成熟的统计推断方法可供我们选用等.总体信息是很重要的信息,为了获取此种信息往往耗资巨大.比如,我国为确认国产轴承寿命分布为韦布尔分布前后花了五年时间,处理了几千个数据后才定下的.

2. 样本信息

样本信息即抽取样本所得观测值提供的信息.譬如,在有了样本观测值后,我们可以根据它大概知道总体的一些特征数,如总体均值、总体方差等在一个什么范围内.这是最“新鲜”的信息,并且越多越好,希望通过样本对总体分布或总体的某些特征作出较精确的统计推断.没有样本就没有统计学可言.

3. 先验信息

如果我们把抽取样本看作做一次试验,则样本信息就是试验中得到的信息.实际中,人们在试验之前对要做的问题在经验上和资料上总是有所了解的,这些信息对统计推断是有益的.先验信息即是抽样(试验)之前有关统计问题的一些信息.一般说来,先验信息来源于经验和历史资料.先验信息在日常生活和工作中是很重要的.先看一个例子.

例 6.5.1

在某工厂的产品中每天要抽检 $n$ 件以确定该厂产品的质量是否满足要求。产品质量可用不合格品率 $p$ 来度量，也可以用 $n$ 件抽查产品中的不合格品件数 $\theta$ 表示。由于生产过程有连续性，可以认为每天的产品质量是有关联的，即是说，在估计现在的 $p$ 时，以前所积累的资料应该是可供使用的，这些积累的历史资料就是先验信息。为了能使用这些先验信息，需要对它进行加工。譬如，在经过一段时间后，就可根据历史资料对过去 $n$ 件产品中的不合格品件数 $\theta$ 构造一个分布

$$P(\theta =i)={\pi }_i,i=1,2,\dots ,n.\text{\hspace{1em}(6.5.1)}$$

这种对先验信息进行加工获得的分布今后称为先验分布.这种先验分布是对该厂过去产品的不合格品率的一个全面看法. 基于上述三种信息进行统计推断的统计学称为贝叶斯统计学.它与经典统计学的差别就在于是否利用先验信息.贝叶斯统计在重视使用总体信息和样本信息的同时，还注意先验信息的收集、挖掘和加工，使它数量化，形成先验分布，参加到统计推断中来.忽视先验信息的利用，有时是一种浪费，有时还会导出不合理的结论. 贝叶斯学派的基本观点是:任一未知量 $\theta$ 都可看作随机变量,可用一个概率分布去描述,这个分布称为先验分布;在获得样本之后,总体分布、样本与先验分布通过贝叶斯公式结合起来得到一个关于未知量 $\theta$ 的新分布——后验分布;任何关于 $\theta$ 的统计推断都应该基于 $\theta$ 的后验分布进行. 关于未知量是否可看作随机变量在经典学派与贝叶斯学派间争论了很长时间.因为任一未知量都有不确定性,而在表述不确定性的程度时,概率与概率分布是最好的语言,因此把它看成随机变量是合理的.如今经典学派已不反对这一观点:著名的美国经典统计学家莱曼(Lehmann,$E$.L.)在他的《点估计理论》一书中写道:“把统计问题中的参数看作随机变量的实现要比看作未知参数更合理一些”.如今两派的争论焦点是:如何利用各种先验信息合理地确定先验分布.这在有些场合是容易解决的,但在很多场合是相当困难的,关于这方面问题的讨论可参阅文献[11].

6.5.2 贝叶斯公式的密度函数形式

贝叶斯公式的事件形式已在 §1.4 节中叙述.这里用随机变量的概率函数再一次叙述贝叶斯公式,并从中介绍贝叶斯学派的一些具体想法. （1）总体依赖于参数 $\theta$ 的概率函数在经典统计中记为 $p(x;\theta )$ ，它表示参数空间 $\Theta$ 中不同的 $\theta$ 对应不同的分布。在贝叶斯统计中应记为 $p(x|\theta )$ ，它表示在随机变量 $\theta$ 取某个给定值时总体的条件概率函数。 (2) 根据参数 $\theta$ 的先验信息确定先验分布 $\pi (\theta )$ . （3）从贝叶斯观点看，样本 $X=(x_1,x_2,\cdots ,x_n)$ 的产生要分两步进行。首先设想从先验分布 $\pi (\theta )$ 产生一个个体 ${\theta }_0$ 。这一步是“老天爷”做的，人们是看不到的，故用“设想”二字。第二步从 $p(X|{\theta }_0)$ 中产生一组样本。这时样本 $X=(x_1,x_2,\cdots ,x_n)$ 的联合条件概率函数为

$$p\left(\text{X}\mid {\theta }_0\right)=p\left(x_1,x_2,\dots ,x_n\mid {\theta }_0\right)=\prod _{i=1}^{n}p\left(x_i\mid {\theta }_0\right),$$

这个分布综合了总体信息和样本信息. （4）由于 ${\theta }_0$ 是设想出来的，仍然是未知的，它是按先验分布 $\pi (\theta )$ 产生的.为把先验信息综合进去，不能只考虑 ${\theta }_0$ ，对 $\theta$ 的其他值发生的可能性也要加以考虑，故要用 $\pi (\theta )$ 进行综合.这样一来，样本 $X$ 和参数 $\theta$ 的联合分布为

$$h(\mathbf{X},\theta )=p(\mathbf{X}\mid \theta )\pi (\theta ).$$

这个联合分布把总体信息、样本信息和先验信息三种可用信息都综合进去了. （5）我们的目的是要对未知参数 $\theta$ 作统计推断。在没有样本信息时，我们只能依据先验分布对 $\theta$ 作出推断。在有了样本观测值 $X=(x_1,x_2,\cdots ,x_n)$ 之后，我们应依据 $h(X,\theta )$ 对 $\theta$ 作出推断。若把 $h(X,\theta )$ 作如下分解：

$$h(\mathbf{X},\theta )=\pi (\theta |\mathbf{X})m(\mathbf{X}),$$

其中 $m(X)$ 是 $X$ 的边际概率函数

$$m(\text{X})={\int }_{\Theta }h(\text{X},\theta )\mathrm{d}\theta ={\int }_{\Theta }p(\text{X}|\mathbf{\theta })\pi (\mathbf{\theta })\mathrm{d}\theta ,\text{\hspace{1em}(6.5.2)}$$

它与 $\theta$ 无关, 或者说 $m(X)$ 中不含 $\theta$ 的任何信息. 因此能用来对 $\theta$ 作出推断的仅是条件分布 $\pi (\theta |X)$ , 它的计算公式是

$$\pi (\theta \mid X)=\frac{h(X,\theta )}{m(X)}=\frac{p(X\mid \theta )\pi (\theta )}{{\int }_{\theta }p(X\mid \theta )\pi (\theta )\mathrm{d}\theta }.\text{\hspace{1em}(6.5.3)}$$

这个条件分布称为 $\theta$ 的后验分布, 它集中了总体、样本和先验中有关 $\theta$ 的一切信息. (6.5.3) 式就是用密度函数表示的贝叶斯公式, 它也是用总体和样本对先验分布 $\pi (\theta )$ 作调整的结果, 它要比 $\pi (\theta )$ 更接近 $\theta$ 的实际情况.

6.5.3 贝叶斯估计

由后验分布 $\pi (\theta |X)$ 估计 $\theta$ 有三种常用的方法： - 使用后验分布的密度函数最大值点作为 $\theta$ 的点估计的最大后验估计. - 使用后验分布的中位数作为 $\theta$ 的点估计的后验中位数估计. - 使用后验分布的均值作为 $\theta$ 的点估计的后验期望估计. 用得最多的是后验期望估计,它一般也简称为贝叶斯估计,记为 ${\hat{\theta }}_B$ .

例 6.5.2

设某事件 $A$ 在一次试验中发生的概率为 $\theta$ ，为估计 $\theta$ ，对试验进行了 $n$ 次独立观测，其中事件 $A$ 发生了 $X$ 次，显然 $X\mid \theta \sim b(n,\theta )$ ，即

$$P\left(X=x\mid \theta \right)=\left(\genfrac{}{}{0ex}{}{n}{x}\right){\theta }^x(1-\theta {)}^{n-x},x=0,1,\dots ,n.$$

假若我们在试验前对事件 $A$ 没有什么了解, 从而对其发生的概率 $\theta$ 也没有任何信息. 在这种场合, 贝叶斯本人建议采用“同等无知”的原则使用区间 $(0,1)$ 上的均匀分布 $U(0,1)$ 作为 $\theta$ 的先验分布, 因为它取 $(0,1)$ 上的每一点的机会均等. 贝叶斯的这个建议被后人称为贝叶斯假设. 由此即可利用贝叶斯公式求出 $\theta$ 的后验分布. 具体如下: 先写出 $X$ 和 $\theta$ 的联合分布

$$h(x,\theta )=\left(\genfrac{}{}{0ex}{}{n}{x}\right){\theta }^x(1-\theta {)}^{n-x},x=0,1,\dots ,n,0<\theta <1,$$

然后求 $X$ 的边际分布

$$m(x)=\left(\genfrac{}{}{0ex}{}{n}{x}\right){\int }_0^1{\theta }^x(1-\theta {)}^{n-x}\mathrm{d}\theta =(\genfrac{}{}{0ex}{}{n}{x})\frac{\Gamma (x+1)\Gamma (n-x+1)}{\Gamma (n+2)},$$

最后求出 $\theta$ 的后验分布

$$\pi (\theta \mid x)=\frac{h(x,\theta )}{m(x)}=\frac{\Gamma (n+2)}{\Gamma (x+1)\Gamma (n-x+1)}{\theta }^{(x+1)-1}(1-\theta {)}^{(n-x+1)-1},0<\theta <1.$$

最后的结果说明 $\theta \mid x\sim Be(x+1,n-x+1)$ ，其后验期望估计为

$${\hat{\theta }}_B=E(\theta \mid x)=\frac{x+1}{n+2}.\text{\hspace{1em}(6.5.4)}$$

假如不用先验信息,只用总体信息与样本信息,那么事件 $A$ 发生的概率的最大似然估计为

$${\hat{\theta }}_M=\frac{x}{n},$$

它与贝叶斯估计是不同的两个估计.某些场合,贝叶斯估计要比最大似然估计更合理一点.比如,在产品抽样检验中只区分合格品和不合格品, $\theta$ 表示不合格品率,对质量好的产品批,抽检的产品常为合格品,但“抽检3个全是合格品”与“抽检10个全是合格品”这两个事件在人们心目中留下的印象是不同的,后者的质量比前者更信得过.这种差别在不合格品率 $\theta$ 的最大似然估计 ${\hat{\theta }}_M$ 中反映不出来(两者都为0),而用贝叶斯估计 ${\hat{\theta }}_B$ 则有所反映,两者分别是 $1/(3+2)=0.20$ 和 $1/(10+2)=0.083$ .类似地,对质量差的产品批,抽检的产品常为不合格品,这时“抽检3个全是不合格品”与“抽检10个全是不合格品”也是有差别的两个事件,前者质量很差,后者则不可救药.这种差别用 ${\hat{\theta }}_M$ 也反映不出(两者都是1),而 ${\hat{\theta }}_B$ 则分别是 $(3+1)/(3+2)=0.80$ 和 $(10+1)/(10+2)=0.917$ .由此可以看到,在这些极端情况下,贝叶斯估计比最大似然估计更符合人们的理念.

例 6.5.3

设 $x_1,x_2,\cdots ,x_n$ 是来自正态分布 $N(\mu ,{\sigma }_0^2)$ 的一个样本，其中 ${\sigma }_0^2$ 已知， $\mu$ 未知，假设 $\mu$ 的先验分布亦为正态分布 $N(\theta ,{\tau }^2)$ ，其中先验均值 $\theta$ 和先验方差 ${\tau }^2$ 均已知，试求 $\mu$ 的贝叶斯估计.

解

样本 $X$ 的分布和 $\mu$ 的先验分布分别为

$$p(\mathbf{X}\mid \mathbf{\mu })=(2\pi {\sigma }_0^2{)}^{-n/2}\exp \left\{-\frac{1}{2{\sigma }_0^2}\sum _{i=1}^n(x_i-\mathbf{\mu }{)}^2\right\},$$$$\pi (\mu )=(2\pi {\tau }^2{)}^{-1/2}\exp \left\{-\frac{1}{2{\tau }^2}(\mu -\theta {)}^2\right\},$$

由此可以写出 $X$ 与 $\mu$ 的联合分布

$$h(X,\mu )=k_1\cdot \exp \left\{-\frac{1}{2}\left[\frac{n{\mu }^2-2n\mu \bar{x}+{\sum }_{i=1}^n{x}_i^2}{{\sigma }_0^2}+\frac{{\mu }^2-2\theta \mu +{\theta }^2}{{\tau }^2}\right]\right\}.$$

其中 $\overline{x}=\frac{1}{n}{\sum }_{i=1}^n{x}_i,k_1=(2\pi {)}^{-(n+1)/2}{\tau }^{-1}{\sigma }_0^{-n}$ . 若记

$$A=\frac{n}{{\sigma }_0^2}+\frac{1}{{\tau }^2},B=\frac{n\overline{x}}{{\sigma }_0^2}+\frac{\theta }{{\tau }^2},C=\frac{{\sum }_{i=1}^n{x}_i^2}{{\sigma }_0^2}+\frac{{\theta }^2}{{\tau }^2},$$

则有

$$h(\mathbf{X},\mu )=k_1\exp \left\{-\frac{1}{2}[A{\mu }^2-2B\mu +C]\right\}=k_1\exp \left\{-\frac{(\mu -B/A{)}^2}{2/A}-\frac{1}{2}(C-B^2/A)\right\}.$$

意到 $A$, $B$, $C$ 均与 $\mu$ 无关, 由此容易算得样本的边际密度函数

$$m(\text{X})={\int }_{-\infty }^{\infty }h(\text{X},\mu )\mathrm{d}\mu =k_1\exp \left\{-\frac{1}{2}(C-B^2/A)\right\}(2\pi /A{)}^{1/2},$$

应用贝叶斯公式即可得到后验分布

$$\pi (\mu \mid X)=\frac{h(X,\mu )}{m(X)}=(2\pi /A{)}^{-1/2}\exp \left\{-\frac{1}{2/A}(\mu -B/A{)}^2\right\}.$$

这说明在样本给定后, $\mu$ 的后验分布为 $N(B/A,1/A)$ ，即

$$\mu \mid X\sim N\left(\frac{n\overline{x}{\sigma }_0^{-2}+\theta {\tau }^{-2}}{n{\sigma }_0^{-2}+{\tau }^{-2}},\frac{1}{n{\sigma }_0^{-2}+{\tau }^{-2}}\right),$$

后验均值即为其贝叶斯估计

$$\hat{\mu }=\frac{n/{\sigma }_0^2}{n/{\sigma }_0^2+1/{\tau }^2}\bar{x}+\frac{1/{\tau }^2}{n/{\sigma }_0^2+1/{\tau }^2}\theta .$$

它是样本均值 $\overline{x}$ 与先验均值 $\theta$ 的加权平均. 当总体方差 ${\sigma }_0^2$ 较小或样本量 $n$ 较大时, 样本均值 $\overline{x}$ 的权重较大; 当先验方差 ${\tau }^2$ 较小时, 先验均值 $\theta$ 的权重较大, 这一综合很符合人们的经验.

6.5.4 共轭先验分布

从贝叶斯公式可以看出,整个贝叶斯统计推断只要先验分布确定后就没有理论上的困难.关于先验分布的确定有多种途径,此处我们介绍一类最常用的先验分布类——共轭先验分布.

定义 6.5.1

设 $\theta$ 是总体分布 $p(x;\theta )$ 中的参数, $\pi (\theta )$ 是其先验分布, 若对任意来自 $p(x;\theta )$ 的样本观测值得到的后验分布 $\pi (\theta |X)$ 与 $\pi (\theta )$ 属于同一个分布族, 则称该分布族是 $\theta$ 的共轭先验分布(族).

例 6.5.4

在例6.5.2 中, 我们知道 $(0,1)$ 上的均匀分布就是贝塔分布的一个特例 $Be(1,1)$ , 其对应的后验分布则是贝塔分布 $Be(x+1,n-x+1)$ . 更一般地, 设 $\theta$ 的先验分布是 $Be(a,b),a>0,b>0,a,b$ , 这说明贝塔分布是伯努利试验中成功概率的共轭先验分布. 类似地,由例6.5.3 可以看出,在方差已知时正态总体均值的共轭先验分布是正态分布.$Be(x+a,n-x+b)$ , 这说明贝塔分布是伯努利试验中成功概率的共轭先验分布. 类似地,由例 6.5.3 可以看出,在方差已知时正态总体均值的共轭先验分布是正态分布.

习题6.5

1. 设一页书上的错别字个数服从泊松分布 $P(\lambda ),\lambda$ 有两个可能取值：1.5 和 1.8，且先验分布为

$$P(\lambda =1.5)=0.45,P(\lambda =1.8)=0.55,$$

现检查了一页,发现有3个错别字,试求 $\lambda$ 的后验分布. 2. 设总体为均匀分布 $U(\theta ,\theta +1),\theta$ 的先验分布是均匀分布 $U(10,16)$ . 现有三个观测值: 11.7, 12.1, 12.0. 求 $\theta$ 的后验分布. 3. 设 $x_1,x_2,\cdots ,x_n$ 是来自几何分布的样本，总体分布列为

$$P(X=k\mid \theta )=\theta (1-\theta {)}^k,k=0,1,2,\dots ,$$

$\theta$ 的先验分布是均匀分布 $U(0,1)$ . (1) 求 $\theta$ 的后验分布； (2) 若 4 次观测值为 4,3,1,6, 求 $\theta$ 的贝叶斯估计. 4. 验证: 泊松分布的均值 $\lambda$ 的共轭先验分布是伽马分布. 5. 验证: 正态总体方差(均值已知)的共轭先验分布是倒伽马分布(称 $X$ 服从倒伽马分布, 如果 1/X 服从伽马分布). 6. 设 $x_1,x_2,\cdots ,x_n$ 是来自如下总体的一个样本

$$p(x\mid \theta )=\frac{2x}{{\theta }^2},0<x<\theta .$$

(1) 若 $\theta$ 的先验分布为均匀分布 $U(0,1)$ ，求 $\theta$ 的后验分布； (2) 若 $\theta$ 的先验分布为 $\pi (\theta )=3{\theta }^2,0<\theta <1$ ，求 $\theta$ 的后验分布. 7. 设 $x_1,x_2,\cdots ,x_n$ 是来自如下总体的一个样本

$$p(x\mid \theta )=\theta x^{\theta -1},0<x<1.$$

若取 $\theta$ 的先验分布为伽马分布，即 $\theta \sim Ga(\alpha ,\lambda )$ ，求 $\theta$ 的后验期望估计. 8. 设 $x_1,x_2,\cdots ,x_n$ 是来自均匀分布 $U(0,\theta )$ 的样本， $\theta$ 的先验分布是帕雷托分布，其密度函数为 $\pi (\theta )=\frac{\beta {\theta }_0^{\beta }}{{\theta }^{\beta +1}},\theta >{\theta }_0$ ，其中 $\beta ,{\theta }_0$ 是两个已知的常数. (1) 验证: 帕雷托分布是 $\theta$ 的共轭先验分布; (2) 求 $\theta$ 的贝叶斯估计. 9. 设指数分布 $Exp(\theta )$ 中未知参数 $\theta$ 的先验分布为伽马分布 $Ga(\alpha ,\lambda )$ ，现从先验信息得知：先验均值为 0.0002，先验标准差为 0.01，试确定先验分布。 10. 设 $x_1,x_2,\cdots ,x_n$ 为来自如下幂级数分布的样本，总体分布密度为

$$p\left(x_1;c,\theta \right)=c{x}_1^{c-1}{\theta }^{-c}I\left\{0⩽x_1⩽\theta \right\}(c>0,\theta >0),$$

证明:(1) 若 $c$ 已知, 则 $\theta$ 的共轭先验分布为帕雷托分布; (2) 若 $\theta$ 已知, 则 $c$ 的共轭先验分布为伽马分布. 11. 某人每天早上在汽车站等公共汽车的时间(单位:min)服从均匀分布 $U(0,\theta )$ ，其中 $\theta$ 未知，假设 $\theta$ 的先验分布为

$$\pi (\theta )=\{\begin{array}{ll}192/{\theta }^4& \theta ⩾4,\\ 0,& \theta <4.\end{array}$$

假如此人在三个早上等车的时间分别为 5,3,8 min, 求 $\theta$ 的后验分布. 12. 从正态总体 $N(\theta ,2^2)$ 中随机抽取容量为 100 的样本, 又设 $\theta$ 的先验分布为正态分布, 证明: 不管先验分布的标准差为多少, 后验分布的标准差一定小于 1/5. 13. 设随机变量 $X$ 服从负二项分布, 其概率分布为

$$f(x\mid p)=\left(\genfrac{}{}{0ex}{}{x-1}{k-1}\right)p^k(1-p{)}^{x-k},x=k,k+1,\dots .$$

证明其成功概率 $p$ 的共轭先验分布族为贝塔分布族. 14. 从一批产品中抽检 100 个, 发现 3 个不合格, 假定该产品不合格率 $\theta$ 的先验分布为贝塔分布 $Be(2,200)$ , 求 $\theta$ 的后验分布.