# 方差分析中的贝叶斯估计问题研究Bayesian Estimation Study in Variance Analysis Model

• 全文下载: PDF(391KB)    PP.508-515   DOI: 10.12677/SA.2017.65057
• 下载量: 136  浏览量: 567

Bayesian estimation is a simple parameter estimation method through the already known data to estimate the unknown parameters of the prior distribution method, which makes the bayesian estimation can make full use of prior information to get a more reasonable estimate results. This paper first introduces the background of the bayesian estimation method and train of thought, and then uses bayesian estimation method to get the concrete expression of the parameter estimation in the analysis of variance model, numerical simulations by programming to solve. At the same time, applying bayesian estimation method is analyzed in the analysis of variance model.

1. 引言

1.1. 研究背景及现状

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

$\pi \left(\theta |x,\eta \right)=\frac{f\left(x|\theta \right)\pi \left(\theta |\eta \right)}{\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)}$

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)$

$N=r\ast n$

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

$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 \mathrm{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(\mu \right)={c}_{1}$ , $\pi \left(\alpha \right)={c}_{2}$ , $\pi \left({\sigma }_{\epsilon }^{2}\right)=\frac{1}{{\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 \mathrm{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}$

$\begin{array}{c}\pi \left({\sigma }_{\epsilon }^{2}|y,\mu ,\alpha \right)=\frac{h\left(y,\mu ,\alpha |{\sigma }_{\epsilon }^{2}\right)}{\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 \mathrm{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 \mathrm{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$

$⇒{\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 }$

$⇒{\alpha }_{\mu }={\left({x}^{\prime }x\right)}^{-1}{x}^{\prime }\left(y-A\alpha \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)$ (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$

$⇒{\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$

$⇒v={\left({A}^{\prime }A+\frac{{\sigma }_{\epsilon }^{2}}{{\sigma }_{\alpha }^{2}}{I}_{r}\right)}^{-1}{A}^{\prime }\left(y-x\mu \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 }_{\epsilon }^{2}}{I}_{r}\right)}^{-1}\right)$ (1-3)

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

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

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

$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 \mathrm{exp}\left\{-\frac{1}{2{\sigma }_{\alpha }^{2}}\left({\alpha }^{\prime }\alpha \right)\right\}$

$\begin{array}{c}\pi \left({\sigma }_{\alpha }^{2}|\alpha \right)=\frac{h\left(\alpha |{\sigma }_{\alpha }^{2}\right)}{\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 \mathrm{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 \mathrm{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. 模拟的结果及分析

4. 实例分析

Table 1. The simulated results of bayesian estimation are obtained from a random generated value with 50 and 100

Table 2. The potential change value of skin under four emotional states (mV)

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

$\begin{array}{c}y=\left(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,\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}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\right)\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$

$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$

5. 总结及展望

 [1] Hsu, J.C. (1992) Simultaneous Inference with Respect to the Best Treatment in Block Designs. Journal of the American Statistical, 6, 107-125. [2] Press, S.J. 贝叶斯统计学, 原理, 模型及应用[M]. 廖文, 陈安贵, 等, 译. 北京: 中国统计出版社, 1992: 187-242. [3] 韦来生. 方差分析模型中参数的经验Bayes估计及其优良性问题[N]. 高校应用数学学报, 1995, 12(2): 163-174. [4] 王松桂, 史建红, 尹素菊, 吴密霞. 线性模型引论[M]. 北京: 科学出版社, 2007: 55-74. [5] 李琪琪, 韦来生. 半参数回归模型中参数的Bayes估计[J]. 中国科学技术大学学报, 2010(9): 881-886. [6] 霍涉云, 张伟平, 韦来生. 一类线性模型参数的Bayes估计及其优良性[J]. 中国科学技术大学学报, 2007(7): 773-776, 802.