1. 引言

1.1. 研究背景及现状

贝叶斯理论，是由英国数学家贝叶斯(Thomas Bayes，1701年~1761年)通过一篇论文“有关随机问题的求解”开始出现的，但开始时它的思想并没有得到重视，直到数学家拉普拉斯(Laplace)用贝叶斯方法重新发现了贝叶斯定理，并导出重要的“相继律”，贝叶斯思想才慢慢地受到人们的重视，形成今天的贝叶斯理论，并广泛应用于数学、工程等领域 [1] 。其中在数学领域里，贝叶斯估计是贝叶斯理论的一个重要组成部分，并广泛的应用于参数估计问题中，其中1992年廖文等 [2] 就写出了书籍《贝叶斯统计学、原理、模型及应用》；1995年，未来生 [3] 写了《方差分析中参数的经验贝叶斯估计及其优良性问题》的论文，介绍经验贝叶斯方法用在了贝叶斯估计中；2007年，王松桂 [4] 等出书《线性模型引论》，介绍线性模型的基础估计研究，并有细致推导；2010年，李琪琪，未来生等 [5] 发表论文《半参数回归模型中参数的贝叶斯估计》，将贝叶斯估计运用到了半参数回归模型。

由于贝叶斯估计是先通过分析已知的过去数据得到先验分布，再通过先验分布得到参数的后验分布，因此，通过贝叶斯估计得到的参数估计值会比一般估计值更精确、更合理，这一理论已经被很多统计学者通过模拟比较得到证实，2007年，霍涉云，张伟平等 [6] 发表论文《一类线性模型参数的贝叶斯估计及其优良性》。因此在本文就不做比较论证。本文主要是将贝叶斯估计运用到方差分析模型的参数估计中，并将理论与Gibbs抽样方法结合通过数值模拟得到相对精确是贝叶斯估计值。同时，将此方法运用到实际问题中进行实例分析。

1.2. 贝叶斯估计的预备知识

贝叶斯估计运用在参数估计中的思路如下：

思路：假设有一组数据 $x=\left(x_1,x_2,\cdots ,x_n\right)$ ，其样本的似然函数为 $f\left(x|\theta \right)$ ，其中参数是 $\theta =\left({\theta }_1,{\theta }_2,\cdots ,{\theta }_k\right)$ ，经典统计学派的思想是：假设 $\theta$ 是可由 $x$ 估计的参数向量；而对于贝叶斯学派的观点是：认为参数向量 $\theta$ 本身就服从先验分布 $\pi \left(\theta |\eta \right)$ 且 $\eta$ 被称为是超参数，得到：

$\pi \left(\theta |x,\eta \right)=\frac{f\left(x|\theta \right)\pi \left(\theta |\eta \right)}{{\displaystyle \int f\left(x|v\right)\pi \left(v|\eta \right)\text{d}v}}=\frac{f\left(x|\theta \right)\pi \left(\theta |\eta \right)}{m\left(x|\eta \right)}$

其中 $m\left(x|\eta \right)$ 是 $x$ 的边际分布， $\pi \left(\theta |x,\eta \right)$ 是后验分布。

这时作为 $\theta$ 的估计可选用后验分布 $\pi \left(\theta |x,\eta \right)$ 的某个位置特征量(即：最大值，中位数和期望值)。其中，由后验分布的任意一个位置特征量求得的估计都被称为是 $\theta$ 的贝叶斯估计，记为 ${\hat{\theta }}_B$ ，有时也被记为 $\hat{\theta }$ 。(注：在一般情况下这三种贝叶斯估计是不同的，所以使用时是按照实际情况选择其中更合适的估计，但是由于正态分布的对称性可以知道，对于正态分布来说， $\theta$ 的三种贝叶斯估计值都是相等的)。

2. 模型系数的贝叶斯估计

将线性模型矩阵化如下：

$\left(\begin{array}{c}{y}_{11}\\ ⋮\\ y_{1n}\\ ⋮\\ y_{r1}\\ ⋮\\ y_{rn}\end{array}\right)=\left(\begin{array}{c}{1}_{n_1}\\ 1_{n_2}\\ ⋮\\ ⋮\\ ⋮\\ 1_{n_{r-1}}\\ 1_{n_r}\end{array}\right)\ast \mu +\left(\begin{array}{ccccccc}{1}_n& 0& \cdots & \cdots & \cdots & 0& 0\\ 0& 1_n& 0& \cdots & \cdots & 0& ⋮\\ 0& 0& \ddots & \ddots & & ⋮& ⋮\\ ⋮& ⋮& \ddots & \ddots & \ddots & ⋮& ⋮\\ ⋮& ⋮& & \ddots & \ddots & 0& ⋮\\ 0& 0& \cdots & \cdots & 0& 1_n& 0\\ 0& 0& \cdots & \cdots & \cdots & 0& 1_n\end{array}\right)\ast \left(\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ ⋮\\ ⋮\\ ⋮\\ {\alpha }_{r-1}\\ {\alpha }_r\end{array}\right)+\left(\begin{array}{c}{\epsilon }_{11}\\ ⋮\\ {\epsilon }_{1n}\\ ⋮\\ {\epsilon }_{r1}\\ ⋮\\ {\epsilon }_{rn}\end{array}\right)$

其中 $1_i={\left(1,1,\cdots ,1\right)}^{\prime }$ 是一个 $i\times 1$ 的全 $1$ 列向量，且 $n_1=n_2=\cdots =n_{r-1}=n_r=r$

设 $N=r\ast n$ 。

故将单因素方差分析模型表示成了线性模型的一般形式：

$y=x\mu +A\alpha +\epsilon$

其中 $\epsilon ~N\left(0,{\sigma }_{\epsilon }^2{I}_N\right)$ ， $\alpha ~N\left(0,{\sigma }_{\alpha }^2{I}_r\right)$ 且 $\mu$ 和 $\alpha =\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_r\right)$ 为参数，则 $y~N\left(x\mu +A\alpha ,{\sigma }_{\epsilon }^2{I}_N\right)$ 。

故 $y$ 在参数 $\mu ,\alpha ,{\sigma }_{\epsilon }^2$ 下的似然函数为：

$L\left(y|\mu ,\alpha ,{\sigma }_{\epsilon }^2\right)=\frac{1}{{\left(\sqrt{2\text{π}}{\sigma }_{\epsilon }\right)}^N}\ast \exp \left\{-\frac{1}{2{\sigma }_{\epsilon }^2}{\left(y-x\mu -A\alpha \right)}^{\prime }\left(y-x\mu -A\alpha \right)\right\}$

假设参数 ${\sigma }_{\epsilon }^2,\mu ,\alpha$ 相互独立，且参数的先验分布函数为：

$\pi \left(\mu \right)=c_1$ , $\pi \left(\alpha \right)=c_2$ , $\pi \left({\sigma }_{\epsilon }^2\right)=\frac{1}{{\sigma }_{\epsilon }^2}$

其中 $c_1$ 和 $c_2$ 为常数。

则参数 ${\sigma }_{\epsilon }^2$ 的联合密度函数为：

$\begin{array}{c}h\left(y,\mu ,\alpha |{\sigma }_{\epsilon }^2\right)=\frac{1}{{\sigma }_{\epsilon }^2}\ast L\left(y|\mu ,\alpha ,{\sigma }_{\epsilon }^2\right)\\ =\frac{1}{{\sigma }_{\epsilon }^2}\ast \frac{1}{{\left(\sqrt{\text{2π}}{\sigma }_{\epsilon }\right)}^N}\ast \exp \left\{-\frac{1}{2{\sigma }_{\epsilon }^2}{\left(y-x\mu -A\alpha \right)}^{\prime }\left(y-x\mu -A\alpha \right)\right\}\end{array}$

由贝叶斯公式可求得参数 ${\sigma }_{\epsilon }^{\text{2}}$ 的联合后验分布为：

$\begin{array}{c}\pi \left({\sigma }_{\epsilon }^2|y,\mu ,\alpha \right)=\frac{h\left(y,\mu ,\alpha |{\sigma }_{\epsilon }^2\right)}{{\displaystyle \int h\left(y,\mu ,\alpha |{\sigma }_{\epsilon }^2\right)\text{d}{\sigma }_{\epsilon }^2}}\\ \propto {\left({\sigma }_{\epsilon }^2\right)}^{-\left(\frac{N}{2}+1\right)}\ast \exp \left\{-\frac{1}{{\sigma }_{\epsilon }^2}\left[\frac{{\left(y-x\mu -A\alpha \right)}^{\prime }\left(y-x\mu -A\alpha \right)}{2}\right]\right\}\end{array}$ (1-1)

(1-1)式恰好为倒Gamma分布密度函数 $\frac{1}{{\sigma }_{\epsilon }^2}\ast \frac{1}{{\left(\sqrt{\text{2π}}{\sigma }_{\epsilon }\right)}^N}\ast \exp \left\{-\frac{1}{2{\sigma }_{\epsilon }^2}{\left(y-x\mu -A\alpha \right)}^{\prime }\left(y-x\mu -A\alpha \right)\right\}$ 的核即：

$\pi \left({\sigma }_{\epsilon }^2|y,\mu ,\alpha \right)~IG\left(\frac{N}{2},\frac{1}{2}{\left(y-x\mu -A\alpha \right)}^{\prime }\left(y-x\mu -A\alpha \right)\right)$

又由正态变量的二次型可知有：

${\mu }^{\prime }{x}^{\prime }{\left({\sigma }_{\epsilon }^2{I}_n\right)}^{-1}x\mu ={\mu }^{\prime }{\Sigma }_{\mu }^{-1}\mu$

$\Rightarrow {\Sigma }_{\mu }={\sigma }_{\epsilon }^2{\left(x^{\prime }x\right)}^{-1}$

将 ${\Sigma }_{\mu }$ 代入下式：

${\mu }^{\prime }{x}^{\prime }{\left({\sigma }_{\epsilon }^2{I}_{n+r}\right)}^{-1}\left(y-A\alpha \right)={\mu }^{\prime }{\Sigma }_{\mu }^{-1}{\alpha }_{\mu }$

$\Rightarrow {\alpha }_{\mu }={\left(x^{\prime }x\right)}^{-1}{x}^{\prime }\left(y-A\alpha \right)$

得到参数 $\mu$ 的满条件后验分布为：

$\pi \left(\mu |y,\alpha ,{\sigma }_{\epsilon }^2\right)~N\left({\left(x^{\prime }x\right)}^{-1}{x}^{\prime }\left(y-A\alpha \right),{\sigma }_{\epsilon }^2{\left(x^{\prime }x\right)}^{-1}\right)$ (1-2)

又由正态变量的二次型可知有：

${\alpha }^{\prime }{A}^{\prime }{\left({\sigma }_{\epsilon }^2{I}_n\right)}^{-1}A\alpha +{\alpha }^{\prime }{\left({\sigma }_{\alpha }^2{I}_r\right)}^{-1}\alpha ={\alpha }^{\prime }{\Sigma }_{\alpha }^{-1}\alpha$

$\Rightarrow {\Sigma }_{\alpha }={\left(\frac{1}{{\sigma }_{\epsilon }^2}{A}^{\prime }A+\frac{1}{{\sigma }_{\alpha }^2}{I}_r\right)}^{-1}$

将 ${\Sigma }_{\alpha }$ 代入下式：

${\alpha }^{\prime }{A}^{\prime }{\left({\sigma }_{\epsilon }^2{I}_{n+r}\right)}^{-1}\left(y-x\mu \right)={\alpha }^{\prime }{\Sigma }_{\alpha }^{-1}v$

$\Rightarrow v={\left(A^{\prime }A+\frac{{\sigma }_{\epsilon }^2}{{\sigma }_{\alpha }^2}{I}_r\right)}^{-1}{A}^{\prime }\left(y-x\mu \right)$

其中 $v=D\alpha$ ，且 $0<D\alpha <1$ 。

得到参数 $\alpha$ 的满条件后验分布为：

$\pi \left(\alpha |y,{\sigma }_{\epsilon }^2,{\sigma }_{\alpha }^2\right)~N\left({\left(A^{\prime }A+\frac{{\sigma }_{\epsilon }^2}{{\sigma }_{\alpha }^2}{I}_r\right)}^{-1}{A}^{\prime }\left(y-x\mu \right),{\left(\frac{1}{{\sigma }_{\epsilon }^2}{A}^{\prime }A+\frac{1}{{\sigma }_{\epsilon }^2}{I}_r\right)}^{-1}\right)$ (1-3)

从上面的结果得到，参数 $$ 和 $\alpha$ 的满条件后验分布都是不常见的概率分布，本文采用Gibbs抽样的方法来抽取。

又 $\alpha$ 在参数 ${\sigma }_{\alpha }^2$ 下的似然函数为：

$L\left(\alpha |{\sigma }_{\alpha }^2\right)=\frac{1}{{\left(\sqrt{2\text{π}}{\sigma }_{\alpha }\right)}^r}\ast \exp \left\{-\frac{1}{2{\sigma }_{\alpha }^2}\left({\alpha }^{\prime }\alpha \right)\right\}$

假设参数 ${\sigma }_{\alpha }^2$ 的先验分布函数为：

$\pi \left({\sigma }_{\alpha }^2\right)=\frac{1}{{\sigma }_{\alpha }^2}$

则参数 $\alpha$ 的联合密度函数为：

$h\left(\alpha |{\sigma }_{\alpha }^2\right)=\frac{1}{{\sigma }_{\alpha }^2}\ast L\left(\alpha |{\sigma }_{\alpha }^2\right)=\frac{1}{{\sigma }_{\alpha }^2}\ast \frac{1}{{\left(\sqrt{\text{2π}}{\sigma }_{\alpha }\right)}^r}\ast \exp \left\{-\frac{1}{2{\sigma }_{\alpha }^2}\left({\alpha }^{\prime }\alpha \right)\right\}$

又由贝叶斯估计的公式可求得参数 ${\sigma }_{\alpha }^{\text{2}}$ 的联合后验分布函数为：

$\begin{array}{c}\pi \left({\sigma }_{\alpha }^2|\alpha \right)=\frac{h\left(\alpha |{\sigma }_{\alpha }^2\right)}{{\displaystyle \int h\left(\alpha |{\sigma }_{\alpha }^2\right)\text{d}{\sigma }_{\alpha }^2}}\\ \propto \frac{1}{{\sigma }_{\alpha }^2}\ast \frac{1}{{\left(\sqrt{\text{2π}}{\sigma }_{\alpha }\right)}^r}\ast \exp \left\{-\frac{1}{2{\sigma }_{\alpha }^2}\left({\alpha }^{\prime }\alpha \right)\right\}\\ \propto {\left({\sigma }_{\alpha }^2\right)}^{-\left(\frac{r}{2}+1\right)}\ast \exp \left\{-\frac{1}{{\sigma }_{\alpha }^2}\left(\frac{{\alpha }^{\prime }\alpha }{2}\right)\right\}\end{array}$ (1-4)

(1-4)式恰好为倒Gamma分布密度函数的核即：

$\pi \left({\sigma }_{\alpha }^2|\alpha \right)~IG\left(\frac{r}{2},\frac{1}{2}{\alpha }^{\prime }\alpha \right)$

故由上可得方差分析模型中参数的后验分布函数为：

$\pi \left({\sigma }_{\alpha }^2|\alpha \right)~IG\left(\frac{r}{2},\frac{1}{2}{\alpha }^{\prime }\alpha \right)$

$\pi \left({\sigma }_{\epsilon }^2|y,\mu ,\alpha \right)~IG\left(\frac{N}{2},\frac{1}{2}{\left(y-x\mu -A\alpha \right)}^{\prime }\left(y-x\mu -A\alpha \right)\right)$

$\pi \left(\mu |y,\alpha ,{\sigma }_{\epsilon }^2\right)~N\left({\left(x^{\prime }x\right)}^{-1}{x}^{\prime }\left(y-A\alpha \right),{\sigma }_{\epsilon }^2{\left(x^{\prime }x\right)}^{-1}\right)$

$\pi \left(\alpha |y,{\sigma }_{\epsilon }^2,{\sigma }_{\alpha }^2\right)~N\left({\left(A^{\prime }A+\frac{{\sigma }_{\epsilon }^2}{{\sigma }_{\alpha }^2}{I}_r\right)}^{-1}{A}^{\prime }\left(y-x\mu \right),{\left(\frac{1}{{\sigma }_{\epsilon }^2}{A}^{\prime }A+\frac{1}{{\sigma }_{\alpha }^2}{I}_r\right)}^{-1}\right)$

3. 方差分析中的贝叶斯估计的数值模拟

3.1. 算法步骤

1) 分别给定参数 $\mu$ 和 $\alpha$ 的初值，以及参数假定 $\epsilon ~N\left(0,{\sigma }_{\epsilon }^2{I}_n\right)$ ；

2) 随机生成一组随机误差的值 $\epsilon ~N\left(0,1\right)$ ；

3) 通过单因素方差分析模型： $y=x\mu +A\alpha +\epsilon$ ，得到1组随机生成的 $y$ 值；

4) 此时，在方差分析模型中我们已随机生成一组样本 $y$ ，下面通过样本求解参数 $\mu ,\alpha$ 和 ${\sigma }_{\epsilon }^2$ 的估计值；

5) 用MCMC方法求解贝叶斯的参数估计问题，其中用Gibbs抽样方法从参数 ${\sigma }_{\alpha }^{\text{2}}$ 和 ${\sigma }_{\alpha }^{\text{2}}$ 的满条件分布核中抽样，不断更新参数 $\mu ,\alpha ,{\sigma }_{\alpha }^{\text{2}}$ 和 ${\sigma }_e^2$ ，重复此步骤 $n$ 次，得到 $n$ 组参数的抽样值，取出后面 $n-m$ 项并分别对其求均值输出，此处 $n$ 和 $m$ 分别取5000和2500；

6) 将(2)~(5)步骤重复 $s$ 次，则可得到 $s$ 组有关参数的抽样均值；

7) 将这 $s$ 组参数的抽样均值再求期望，并与原假设数据做均方差矩阵MSEM的值；

8) 分别取 $s$ 等于10和50求解参数估计的值。

3.2. 模拟的结果及分析

设定单因素方差分析模型： $y=x\mu +A\alpha +\epsilon$ 中的参数的真值分别为： $\alpha ={\left({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4\right)}^{\prime }={\left(1,0.2,2,0.6,0,0.8,0.1,1,2,0,-2,-0.6,-0.1,0,-0.02,-0.8,-1,-1,0,-2\right)}^{\prime }$ ， $\mu =5$ 设 $x_{\left(i\right)}$ 是 $i\ast 1$ 的全1列向量，则 $x$ 为 $x_{\left(2000\right)}$ 即为 $2000\ast 1$ 的全1列向量， $A$ 是 $2000\ast 20$ 的矩阵有A = kronecker(diag(20)， $x_{\left(2000\right)}$ )，再随机生成1个 $\epsilon =\left({\epsilon }_1,{\epsilon }_2,\cdots ,{\epsilon }_{2000}\right)$ 其中有 ${\epsilon }_i~N\left(0,1\right)$ ，循环上面过程 $s$ 次，得到 $s$ 组 $y$ 值，再用每组 $y$ 值经过5000次抽样求解参数 $\mu ,\alpha$ 和 ${\sigma }_e$ 的贝叶斯估计，最后求得估计均值见表1。

表1充分的显示出贝叶斯估计与真值的近似程度很高，并且具有很好的优良性和有效性。

4. 实例分析

为研究人们在催眠状态下对各种情绪的反应是否有差异，选取了8个受试者。在催眠状态下，要求每人按任意次序做出恐惧，愉快，忧虑和平静4种反应。表2给出了各受试者在处于这4种情绪状态下皮肤的电位变化值。试在 $\alpha =0.05$ 下，检验受试者在催眠状态下对这4种情绪的反应力是否有显著差异。

表1. 对50和100个随机生成的 $y$ 值求解贝叶斯估计的模拟结果

表2. 4种情绪状态下皮肤的电位变化值(mV)

对于上面的问题建立方差分析模型：

$y=x\mu +A\alpha +\epsilon$

其中：

$\begin{array}{c}y=(23.1,57.6,10.5,23.6,11.9,54.6,21.0,20.3,22.7,53.2,9.7,19.6,13.8,47.1,13.6,23.6,\\ 22.5,53.7,10.8,21.1,13.7,39.2,13.7,16.3,22.6,53.1,8.3,21.6,13.3,37.0,14.8,14.8)\end{array}$

$x$ 为32 × 1的全1列向量，A为kronecker(diag(4),array(d,dim=c(8,1)))一个32 × 4的矩阵，则下面通过R程序得以下结果：

$\mu =25.1489461$ , ${\sigma }_{\epsilon }^2=15.7876821$

${\alpha }_1=0.05077693$ , ${\alpha }_2=-0.011519118$ , ${\alpha }_3=-0.005867851$ , ${\alpha }_4=-0.008307687$

下面检验受试者在催眠状态下对这4种情绪的反应力是否有显著差异， $\alpha =0.05$

假设原假设为： $H_0:{\alpha }_1-{\alpha }_2={\alpha }_1-{\alpha }_3={\alpha }_1-{\alpha }_4=0$

则可知：

$p\left({\alpha }_1\text{,}{\alpha }_2\right)=0.01513101<0.025$

$p\left({\alpha }_1\text{,}{\alpha }_3\right)=0.04791839<0.025$

$p\left({\alpha }_1\text{,}{\alpha }_4\right)=0.06331945>0.025$

由上式可知在 $\alpha =0.05$ 的拒绝区间下，应拒绝原假设，得到结论是受试者在催眠状态下对这4种情绪的反应力有显著差异。

5. 总结及展望

通过将贝叶斯估计运用到单因素方差分析模型中的数值模拟，可以发现贝叶斯估计的优良性和高度的近似性，用在实例分析中也简单易懂。再通过对单因素方差分析模型的研究，可以将贝叶斯估计推广到多因素方差分析模型的研究中。同时也可以经验贝叶斯估计运用在单因素方差分析和多因素方差分析模型的研究，也将是一个新的挑战。

参考文献