CN102521832A - 图像分析方法和系统 - Google Patents

Abstract

一种图像分析方法，包括如下步骤：获取一组T1加权图像，包括对照组图像和实验组图像；对所有T1加权图像预处理；将预处理后的T1加权图像非线性配准到预设的模板图像中；将非线性配准后的T1加权图像进行分割，产生多个感兴趣区域；提取每个感兴趣区域的至少一种特征数据，得到特征矩阵信息；根据预设的贝叶斯多层模型和所述特征矩阵信息，计算得到一参数值；根据所述参数值与预设阀值的关系，判断实验组图像和对照组图像之间感兴趣区域的差异显著性。本发明还提供一种对应的图像分析系统。上述图像分析方法和系统相对传统技术可提取的待分析参数更多更丰富，信息更多样化，能同时处理多区域、多特征的显著性统计检验，使形态学分析更全面。

Description

图像分析方法和系统

【技术领域】

本发明涉及磁共振计算解剖学领域，特别是涉及一种对磁共振图像的分析方法和系统。

【背景技术】

人脑的发育几乎就是人类几千万年大脑神经组织发展进化的重演，它也有着一个由低到高的生理发育过程和生物发展规律。人脑的发展可以划分为三个时期：从出生到成年的脑的加速成长期；脑发育成熟后的相对稳定期；人到老年后的脑的衰老期。由于人脑结构在自身的发展过程中形成了一定的规律，人脑成熟或衰退的部位、时间和表现也不一样。认识大脑的发展规律对认识人类的智力、生理、情感发展有很大的帮助；同时，认识了大脑的发展规律也能帮助人类了解更多的脑部病理变化，如老年痴呆，癫痫，抑郁症等等。无论是大脑的发育衰老，还是脑部的神经性病变，都会引起大脑的解剖形态变化。研究人员和临床医师均希望能绘制出脑部的形变图谱，以帮助了解大脑结构形态变化的区域和过程。形变图谱能反映大脑从发育成熟到衰老的结构纵向变化过程，能反映某一疾病的患者脑部与正常人脑部结构的平均差异，能反映疾病患者随时间发展脑部的病理改变情况，还能反映治疗过程中患者在一段时间内形变的控制程度。同时，大脑各解剖区域与人体生理功能和心理状态存在对应关系，因此，形变图谱所反映的脑部结构形态信息，能帮助研究人员或临床医师获得或预测人体相应的生理、心理状态情况。

在磁共振计算解剖学领域，研究人员常用统计分析法对大脑解剖结构进行形态学分析，确定大脑的形变部位，帮助了解脑部的病变情况或治疗效果。当前常用的形态学分析法包括基于体素的分析法(VBM，Voxel-BasedMorphometry)、基于形变的分析法(DBM，Deformation-Based Morphometry)，以及基于张量的分析法(TBM，Tensor-Based Morphometry)。其中DBM，TBM利用了图像配准得到的形变场进行统计分析，VBM则利用了配准残差。

几种形态学分析法对图像的前期处理采取了相同的做法。首先，选取一幅模板图像作为标准空间，该模板图像或是领域研究人员公认的开放模板(如MNI，Talairach等)，或是根据手头数据独自建立的特定模板；然后将所有待处理图像进行灰度校正，去头骨等处理，并将预处理后的结果线性配准到模板图像空间。在提取统计分析对象时，VBM与DBM、TBM方法各异。在将预处理结果线性配准到模板空间后，VBM将大脑分割成子区域，将目标子区域平滑并规范化，平滑规范化后的结果作为统计分析的输入；而DBM、TBM在得到预处理结果后继续将结果非线性配准到模板图像空间，得到每幅图像对应的配准形变场，对形变场完成一定的变换后作为统计分析的输入。在统计分析时，t检验(t-test)和置换检验(permutation test)是常用的手段。在分析时，每个体素被看成一个独立的检验，因此目标区域体素的个数便是独立检验的次数。

在现有方法进行统计分析时，尽管前期处理对图像或形变场进行了平滑，但对每个像素进行独立的统计检验，仍然丢弃了大脑子区域内的关联性。我们知道大脑有功能分区，每个小的功能区域有相应的生理作用，在病理改变时往往一个区域具有相同的变化趋势，将这样一个子区域作为整体分析更容易对结果进行合理的解释。而前期的平滑处理也只是体素空间位置上的信息叠加，而不是潜在功能区域的信息融合。

此外，从统计学角度看，t检验、置换检验等一般统计学检验方法不太适合小样本大规模回归量的分析，也就是说，由于可用于研究的图像数量有限，对于一个体素或一个子区域，用于描述该目标的特征的维度不能太高。这在一定程度上限制了形态学分析的范围，现有方法只对形态特征描述的其中某一项进行分析，如大脑皮质厚度，沟回宽度等，从而导致分析结果不全面。

【发明内容】

基于此，有必要提供一种能使分析结果更全面的图像分析方法和系统。

一种图像分析方法，包括如下步骤：

获取一组T1加权图像，包括对照组图像和实验组图像；

对所有T1加权图像预处理；

将预处理后的T1加权图像非线性配准到预设的模板图像中；

将非线性配准后的T1加权图像进行分割，产生多个感兴趣区域；

提取每个感兴趣区域的至少一种特征数据，得到特征矩阵信息；

根据预设的贝叶斯多层模型和所述特征矩阵信息，计算得到一参数值，所述预设的贝叶斯多层模型如下：

Pr{Y＝1|C，γ(X)}＝Φ{α^TC+γ(X)}，

γ(X)～MN{0，τ^-1K(X)}，

其中，Y_i为二进制响应变量，Y_i＝1表示第i幅图像来自实验组，Y_i＝0表示第i幅图像来自对照组，i∈{1，2，…，n}，Y_n×1为二进制响应结果，C_n×1为协变量，X_n×p为所述特征矩阵信息，Φ(.)为标准正态累积分布函数，α为回归系数，γ(X)＝{γ(x₁)，γ(x₂)，…，γ(x_n)}表示均值为0，协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)的高斯随机过程，τ服从Gamma分布，K是一个n×n矩阵，K_ij＝K(x_i，x_j)，K为核函数，所述参数值为所述预设的贝叶斯多层模型中τ的值；

根据所述参数值与预设阀值的关系，判断实验组图像和对照组图像之间感兴趣区域的差异显著性。

一种图像分析系统，包括获取单元、处理单元、配准单元、分割单元、提取单元、计算单元和分析单元：

获取单元用于获取一组T1加权图像，包括对照组图像和实验组图像；

处理单元用于对所有T1加权图像预处理；

配准单元用于将预处理后的T1加权图像非线性配准到预设的模板图像中；

分割单元用于将非线性配准后的T1加权图像进行分割，产生多个感兴趣区域；

提取单元用于提取每个感兴趣区域的至少一种特征数据，得到特征矩阵信息；

计算单元用于根据预设的贝叶斯多层模型和所述特征矩阵信息，计算得到一参数值，所述预设的贝叶斯多层模型如下：

Pr{Y＝1|C，γ(X)}＝Φ{α^TC+γ(X)}，

γ(X)～MN{0，τ^-1K(X)}，

其中，Y_i为二进制响应变量，Y_i＝1表示第i幅图像来自实验组，Y_i＝0表示第i幅图像来自对照组，i∈{1，2，…，n}，Y_n×1为二进制响应结果，C_n×1为协变量，X_n×p为所述特征矩阵信息，Φ(.)为标准正态累积分布函数，α为回归系数，γ(X)＝{γ(x₁)，γ(x₂)，…，γ(x_n)}表示均值为0，协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)的高斯随机过程，τ服从Gamma分布，K是一个n×n矩阵，K_ij＝K(x_i，x_j)，K为核函数，所述参数值为所述预设的贝叶斯多层模型中τ的值；

分析单元用于根据所述参数值与预设阀值的关系，判断实验组图像和对照组图像之间感兴趣区域的差异显著性

本发明一较佳实施例中，所述提取每个感兴趣区域的至少一种特征数据，得到特征矩阵信息步骤是从感兴趣区域R_q，q∈{1，2，…，Q}中分别提取每幅图像的多种特征信息：

其中i为图像索引号i∈{1，2，…，n}，q为感兴趣区域索引号q∈{1，2，…，Q}，p为该区域特征信息的维度p≥1，对每个感兴趣区域R_q得到的所述特征矩阵信息如下：

其中x_ij表示第i幅图像感兴趣区域的第j个参数元素，i∈{1，2，…，n}，j∈{1，2，…，p}。

本发明一较佳实施例中，所述预设的贝叶斯多层模型计算方式如下：

首先，定义了一个隐式变量Z＝(z₁，z₂，…，z_n)，Z～MN{α^TC+γ(X)，I}，使得

$Y_i=\left\{\begin{array}{cc}1& \mathrm{if}{Z}_i\&GreaterEqual;0\\ 0& \mathrm{if}{Z}_i<0\end{array}\right.,$

将协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)中的K(x，x′)定义为高斯核，则：

$K(x,x^{\&prime;})=e^{-\frac{\left|\right|x-x^{\&prime;}\left|\right|}{2}},$

然后，假设γ(X)～MN{0，τ^-1K(X)}，α～MN(0，φ^-1)和τ～Gamma(a，b)，则联合后验分布为：

$[Z,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;},\mathrm{\&tau;}|y,C,X]\&Proportional;\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}(I_{Z_i\&GreaterEqual;0}{I}_{y_i=1}+I_{Z_i<0}{I}_{y_i=0})\&times;\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1)$

$\&times;\mathrm{MN}\{\mathrm{\&gamma;};0,{\mathrm{\&tau;}}^{-1}K\left(X\right)\}\&times;\mathrm{Gamma}(\mathrm{\&tau;};a,b)$

$\&times;N({\mathrm{\&alpha;}}_0;0,{\mathrm{\&phi;}}^{-1})\&times;N({\mathrm{\&alpha;}}_1;0,{\mathrm{\&phi;}}^{-1}),$

得到Z_i的满条件概率分布的闭合形式：

$\left[Z_i\right|y_i=1,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,{\mathrm{\&gamma;}}_i,\mathrm{\&tau;}]\&Proportional;I_{Z_i<0}{I}_{y_i=0}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1)$

$\left[Z_i\right|y_i=0,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,{\mathrm{\&gamma;}}_i,\mathrm{\&tau;}]\&Proportional;I_{Z_i\&GreaterEqual;0}{I}_{y_i=1}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1),$

该分布为截断正态分布。α，γ，τ的满条件分布求得如下：

$\left[{\mathrm{\&alpha;}}_0\right|Z_i,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;},\mathrm{\&tau;}]\&Proportional;N\{\frac{{\mathrm{\&Sigma;}}_{i=1}^n(Z_i-{\mathrm{\&alpha;}}_1{C}_i-{\mathrm{\&gamma;}}_i)}{n+\mathrm{\&phi;}},{(n+\mathrm{\&phi;})}^{-1}\}$

$\left[{\mathrm{\&alpha;}}_1\right|Z_i,{\mathrm{\&alpha;}}_0,\mathrm{\&gamma;},\mathrm{\&tau;}]\&Proportional;N\{\frac{{\mathrm{\&Sigma;}}_{i=1}^n{C}_i(Z_i-{\mathrm{\&alpha;}}_0-{\mathrm{\&gamma;}}_i)}{{\mathrm{\&Sigma;}}_{i=1}^n{C}_i^2+\mathrm{\&phi;}},{(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{C}_i^2+\mathrm{\&phi;})}^{-1}\}$

[γ|Z_i，α₀，α₁，τ]∝MN{[1+(τ^-1K(X))^-1](Z-α′C)，[1+τ^-1K(X))^-1]^-1}

$\left[\mathrm{\&tau;}\right|Z_i,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;}]\&Proportional;\mathrm{Gamma}\{a+\frac{1}{2},b+\frac{1}{2}\mathrm{\&gamma;}{\left(K\left(X\right)\right)}^{-1}\mathrm{\&gamma;}\};$

使用Gibbs抽样算法循环计算各个参数的后验分布，Gibbs抽样算法如下：

步骤1：初始化参数[Z⁽⁰⁾，α⁽⁰⁾，γ⁽⁰⁾，τ⁽⁰⁾]；

步骤2：在第t次迭代时，

(i)求Z^(t)：[Z|Z^(t-1)，α^(t-1)，γ^(t-1)]；

(ii)求α^(t)：[α|Z^(t)，α^(t-1)，γ^(t-1)]；

(iii)求γ^(t)：[γ|Z^(t)，α^(t)，τ^(t-1)]；

(iv)求τ^(t)：[τ|γ^(t)]。

步骤3：增加迭代次数t直至到达规定的迭代次数；

通过上述Gibbs抽样算法，得到马尔可夫链蒙特卡洛(MCMC)样本用于贝叶斯推断和预测，计算95％的贝叶斯置信区间并得到点估计量，由此得到τ^-1。

通过上述图像分析方法和系统，把大脑感兴趣区域作为统一整体进行形态学统计，分割的子区域具有一定的功能或解剖结构联系，有助于检测更细微的变化，也有助于提高组织形变与生理行为的对应。与一般的统计检验方式相比，上述图像分析方法和系统对待分析参数的维度没有要求，尤其适用于小样本大规模回归量的分析，可提取的待分析参数更多更丰富，信息更多样化，能同时处理多区域、多特征的显著性统计检验，使形态学分析更全面。。

【附图说明】

图1为本发明一实施例的图像分析方法的步骤流程图；

图2为本发明一实施例的图像分析系统的功能模块图。

【具体实施方式】

为了解决传统图像分析技术分析不全面的问题，提出了一种能使分析结果更全面的图像分析方法和系统。

如图1所示，其为本发明一实施例的图像分析方法的步骤流程图，包括如下步骤：

步骤S101，获取一组T1加权图像，包括对照组图像和实验组图像。

如，获取的T1加权图像为：

其中

表示对照组图像，r∈{1，2，…，L}，表示实验组图像，s∈{1，2，…，M}，n＝L+M。

步骤S102，对所有T1加权图像预处理。预处理根据需求而定，如：灰度校正，去头骨等。

步骤S103，将预处理后的T1加权图像非线性配准到预设的模板图像中。

步骤S104，将非线性配准后的T1加权图像进行分割，产生多个感兴趣区域。如：分割成为Q个感兴趣区域R_q，q∈{1，2，…，Q}。分割动作可以是根据用户输入而分割，也可以是根据预设规则自动分割。

步骤S105，提取每个感兴趣区域的至少一种特征数据，得到特征矩阵信息。

本发明一实施例中，从感兴趣区域R_q，q∈{1，2，…，Q}中分别提取每幅图像的多种特征信息：

其中i为图像索引号i∈{1，2，…，n}，q为感兴趣区域索引号q∈{1，2，…，Q}，p为该区域特征信息的维度p≥1。因此，对每个感兴趣区域R_q，可得到一个数据矩阵如下：

其中x_ij表示第i幅图像感兴趣区域的第j个参数元素，i∈{1，2，…，n}，j∈{1，2，…，p}。

步骤S106，根据预设的贝叶斯多层模型和所述特征矩阵信息，计算得到一参数值，所述预设的贝叶斯多层模型如下：

Pr{Y＝1|C，γ(X)}＝Φ{α^TC+γ(X)}，

γ(X)～MN{0，τ^-1K(X)}，

其中，Y_i为二进制响应变量，Y_i＝1表示第i幅图像来自实验组，Y_i＝0表示第i幅图像来自对照组，i∈{1，2，…，n}。Y_n×1为二进制响应结果，C_n×1为协变量(可描述实验个体的临床参数，如年龄、性别等)，X_n×p为所述特征矩阵信息，Φ(.)为标准正态累积分布函数，α为回归系数，γ(X)＝{γ(x₁)，γ(x₂)，…，γ(x_n)}表示均值为0，协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)的高斯随机过程，τ服从Gamma分布，K是一个n×n矩阵，K_ij＝K(x_i，x_j)，所述参数值为所述预设的贝叶斯多层模型中τ的值。

步骤S107，根据所述参数值与预设阀值的关系，判断实验组图像和对照组图像之间感兴趣区域的差异显著性。

通过上述贝叶斯多层模型，将实验组与对照组感兴趣区域的差异显著性判断转换成检验γ(X)在0处是否为点概率分布，γ(X)是否有一个常数协方矩阵。即，检验τ^-1是否为0，若τ^-1＝0，接受H0，表明实验组与对照组的该感兴趣区域无显著差异，若τ^-1的置信区间远离0，则表明实验组与对照组有显著差异；此外，检验τ^-1K(X)是否等于τ^-1I，若τ^-1K(X)＝τ^-1I，同样接受H₀，表明实验组与对照组无显著差异。

通过上述图像分析方法，把大脑感兴趣区域作为统一整体进行形态学统计，分割的子区域具有一定的功能或解剖结构联系，有助于检测更细微的变化，也有助于提高组织形变与生理行为的对应。与一般的统计检验方法相比，上述图像分析方法对待分析参数的维度没有要求，尤其适用于小样本大规模回归量的分析，可提取的待分析参数更多更丰富，信息更多样化，能同时处理多区域、多特征的显著性统计检验，使形态学分析更全面。

为了进一步提高图像分析效率，本发明还提供一种高效计算上述贝叶斯多层模型的方法，包括如下步骤：

首先，定义了一个隐式变量Z＝(z₁，z₂，…，z_n)，Z～MN{α^TC+γ(X)，I}，使得

$Y_i=\left\{\begin{array}{cc}1& \mathrm{if}{Z}_i\&GreaterEqual;0\\ 0& \mathrm{if}{Z}_i<0\end{array}\right..$

将协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)中的K(x，x′)定义为高斯核，则：

$K(x,x^{\&prime;})=e^{-\frac{\left|\right|x-x^{\&prime;}\left|\right|}{2}}.$

然后，假设γ(X)～MN{0，τ^-1K(X)}，α～MN(0，φ^-1)和τ～Gamma(a，b)，则联合后验分布为：

$[Z,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;},\mathrm{\&tau;}|y,C,X]\&Proportional;\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}(I_{Z_i\&GreaterEqual;0}{I}_{y_i=1}+I_{Z_i<0}{I}_{y_i=0})\&times;\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1)$

$\&times;\mathrm{MN}\{\mathrm{\&gamma;};0,{\mathrm{\&tau;}}^{-1}K\left(X\right)\}\&times;\mathrm{Gamma}(\mathrm{\&tau;};a,b).$

$\&times;N({\mathrm{\&alpha;}}_0;0,{\mathrm{\&phi;}}^{-1})\&times;N({\mathrm{\&alpha;}}_1;0,{\mathrm{\&phi;}}^{-1})$

由于本模型是一个probit回归模型，可得到Z_i的满条件概率分布的闭合形式：

$\left[Z_i\right|y_i=1,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,{\mathrm{\&gamma;}}_i,\mathrm{\&tau;}]\&Proportional;I_{Z_i<0}{I}_{y_i=0}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1),$

$\left[Z_i\right|y_i=0,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,{\mathrm{\&gamma;}}_i,\mathrm{\&tau;}]\&Proportional;I_{Z_i\&GreaterEqual;0}{I}_{y_i=1}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1)$

该分布为截断正态分布。α，γ，τ的满条件分布求得如下：

$\left[{\mathrm{\&alpha;}}_0\right|Z_i,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;},\mathrm{\&tau;}]\&Proportional;N\{\frac{{\mathrm{\&Sigma;}}_{i=1}^n(Z_i-{\mathrm{\&alpha;}}_1{C}_i-{\mathrm{\&gamma;}}_i)}{n+\mathrm{\&phi;}},{(n+\mathrm{\&phi;})}^{-1}\}$

$\left[{\mathrm{\&alpha;}}_1\right|Z_i,{\mathrm{\&alpha;}}_0,\mathrm{\&gamma;},\mathrm{\&tau;}]\&Proportional;N\{\frac{{\mathrm{\&Sigma;}}_{i=1}^n{C}_i(Z_i-{\mathrm{\&alpha;}}_0-{\mathrm{\&gamma;}}_i)}{{\mathrm{\&Sigma;}}_{i=1}^n{C}_i^2+\mathrm{\&phi;}},{(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{C}_i^2+\mathrm{\&phi;})}^{-1}\}.$

[γ|Z_i，α₀，α₁，τ]∝MN{[1+(τ^-1K(X))^-1](Z-α′C)，[1+τ^-1K(X))^-1]^-1}

$\left[\mathrm{\&tau;}\right|Z_i,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;}]\&Proportional;\mathrm{Gamma}\{a+\frac{1}{2},b+\frac{1}{2}\mathrm{\&gamma;}{\left(K\left(X\right)\right)}^{-1}\mathrm{\&gamma;}\}$

由于所有参数的满条件分布均可得到闭合形式，因此使用Gibbs抽样算法循环计算各个参数的后验分布。Gibbs抽样算法如下：

步骤1：初始化参数[Z⁽⁰⁾，α⁽⁰⁾，γ⁽⁰⁾，τ⁽⁰⁾]。

步骤2：在第t次迭代时，

(i)求Z^(t)：[Z|Z^(t-1)，α^(t-1)，γ^(t-1)]；

(ii)求α^(t)：[α|Z^(t)，α^(t-1)，γ^(t-1)]；

(iii)求γ^(t)：[γ|Z^(t)，α^(t)，τ^(t-1)]；

(iv)求τ^(t)：[τ|γ^(t)]。

步骤3：增加迭代次数t直至到达规定的迭代次数。

通过上述Gibbs抽样算法，可得到马尔可夫链蒙特卡洛(MCMC)样本用于贝叶斯推断和预测，计算95％的贝叶斯置信区间，并得到点估计量，由此得到τ^-1，检验τ^-1是否为0或τ^-1K(X)是否等于τ^-1I。

如图2所示，其为本发明一实施例的图像分析系统20的功能模块图，包括：获取单元201、处理单元202、配准单元203、分割单元204、提取单元205、计算单元206和分析单元207。

获取单元201用于获取一组T1加权图像，包括对照组图像和实验组图像。

如，获取的T1加权图像为：

其中

表示对照组图像，r∈{1，2，…，L}，表示实验组图像，s∈{1，2，…，M}，n＝L+M。

处理单元202用于对所有T1加权图像预处理。预处理根据需求而定，如：灰度校正，去头骨等。

配准单元203用于将预处理后的T1加权图像非线性配准到预设的模板图像中。

分割单元204用于将非线性配准后的T1加权图像进行分割，产生多个感兴趣区域。如：分割成为Q个感兴趣区域R_q，q∈{1，2，…，Q}。分割动作可以是根据用户输入而分割，也可以是根据预设规则自动分割。

提取单元205用于提取每个感兴趣区域的至少一种特征数据，得到特征矩阵信息。

本发明一实施例中，从感兴趣区域R_q，q∈{1，2，…，Q}中分别提取每幅图像的多种特征信息：

其中i为图像索引号i∈{1，2，…，n}，q为感兴趣区域索引号q∈{1，2，…，Q}，p为该区域特征信息的维度p≥1。因此，对每个感兴趣区域R_q，可得到一个数据矩阵如下：

其中x_ij表示第i幅图像感兴趣区域的第j个参数元素，i∈{1，2，…，n}，j∈{1，2，…，p}。

计算单元206用于根据预设的贝叶斯多层模型和所述特征矩阵信息，计算得到一参数值，所述预设的贝叶斯多层模型如下：

Pr{Y＝1|C，γ(X)}＝Φ{α^TC+γ(X)}，

γ(X)～MN{0，τ^-1K(X)}，

其中，Y_i为二进制响应变量，Y_i＝1表示第i幅图像来自实验组，Y_i＝0表示第i幅图像来自对照组，i∈{1，2，…，N}。Y_n×1为二进制响应结果，C_n×1为协变量(可描述实验个体的临床参数，如年龄、性别等)，X_n×p为所述特征矩阵信息，Φ(.)为标准正态累积分布函数，α为回归系数，γ(X)＝{γ(x₁)，γ(x₂)，…，γ(x_n)}表示均值为0，协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)的高斯随机过程，τ服从Gamma分布，K是一个n×n矩阵，K_ij＝K(x_i，x_j)，所述参数值为所述预设的贝叶斯多层模型中τ的值。

分析单元207用于根据所述参数值与预设阀值的关系，判断实验组图像和对照组图像之间感兴趣区域的差异显著性。

通过上述贝叶斯多层模型，将实验组与对照组感兴趣区域的差异显著性判断转换成检验γ(X)在0处是否为点概率分布，γ(X)是否有一个常数协方矩阵。即，检验τ^-1是否为0，若τ^-1＝0，接受H₀，表明实验组与对照组的该感兴趣区域无显著差异，若τ^-1的置信区间远离0，则表明实验组与对照组有显著差异；此外，检验τ^-1K(X)是否等于τ^-1I，若τ^-1K(X)＝τ^-1I，同样接受H₀，表明实验组与对照组无显著差异。

通过上述图像分析系统，把大脑感兴趣区域作为统一整体进行形态学统计，分割的子区域具有一定的功能或解剖结构联系，有助于检测更细微的变化，也有助于提高组织形变与生理行为的对应。与一般的统计检验方式相比，上述图像分析系统对待分析参数的维度没有要求，尤其适用于小样本大规模回归量的分析，可提取的待分析参数更多更丰富，信息更多样化，能同时处理多区域、多特征的显著性统计检验，使形态学分析更全面。

为了进一步提高图像分析效率，本发明还提供一种高效计算上述贝叶斯多层模型的方法，包括如下步骤：

首先，定义了一个隐式变量Z＝(z₁，z₂，…，z_n)，Z～MN{α^TC+γ(X)，I}，使得

$Y_i=\left\{\begin{array}{cc}1& \mathrm{if}{Z}_i\&GreaterEqual;0\\ 0& \mathrm{if}{Z}_i<0\end{array}\right..$

将协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)中的K(x，x′)定义为高斯核，则：

$K(x,x^{\&prime;})=e^{-\frac{\left|\right|x-x^{\&prime;}\left|\right|}{2}}.$

然后，假设γ(X)～MN{0，τ^-1K(X)}，α～MN(0，φ^-1)和τ～Gamma(a，b)，则联合后验分布为：

$[Z,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;},\mathrm{\&tau;}|y,C,X]\&Proportional;\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}(I_{Z_i\&GreaterEqual;0}{I}_{y_i=1}+I_{Z_i<0}{I}_{y_i=0})\&times;\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1)$

$\&times;\mathrm{MN}\{\mathrm{\&gamma;};0,{\mathrm{\&tau;}}^{-1}K\left(X\right)\}\&times;\mathrm{Gamma}(\mathrm{\&tau;};a,b).$

$\&times;N({\mathrm{\&alpha;}}_0;0,{\mathrm{\&phi;}}^{-1})\&times;N({\mathrm{\&alpha;}}_1;0,{\mathrm{\&phi;}}^{-1})$

由于本模型是一个probit回归模型，可得到Z_i的满条件概率分布的闭合形式：

$\left[Z_i\right|y_i=1,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,{\mathrm{\&gamma;}}_i,\mathrm{\&tau;}]\&Proportional;I_{Z_i<0}{I}_{y_i=0}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1),$

$\left[Z_i\right|y_i=0,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,{\mathrm{\&gamma;}}_i,\mathrm{\&tau;}]\&Proportional;I_{Z_i\&GreaterEqual;0}{I}_{y_i=1}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1)$

该分布为截断正态分布。α，γ，τ的满条件分布求得如下：

$\left[{\mathrm{\&alpha;}}_0\right|Z_i,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;},\mathrm{\&tau;}]\&Proportional;N\{\frac{{\mathrm{\&Sigma;}}_{i=1}^n(Z_i-{\mathrm{\&alpha;}}_1{C}_i-{\mathrm{\&gamma;}}_i)}{n+\mathrm{\&phi;}},{(n+\mathrm{\&phi;})}^{-1}\}$

$\left[{\mathrm{\&alpha;}}_1\right|Z_i,{\mathrm{\&alpha;}}_0,\mathrm{\&gamma;},\mathrm{\&tau;}]\&Proportional;N\{\frac{{\mathrm{\&Sigma;}}_{i=1}^n{C}_i(Z_i-{\mathrm{\&alpha;}}_0-{\mathrm{\&gamma;}}_i)}{{\mathrm{\&Sigma;}}_{i=1}^n{C}_i^2+\mathrm{\&phi;}},{(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{C}_i^2+\mathrm{\&phi;})}^{-1}\}.$

[γ|Z_i，α₀，α₁，τ]∝MN{[1+(τ^-1K(X))^-1](Z-α′C)，[1+τ^-1K(X))^-1]^-1}

$\left[\mathrm{\&tau;}\right|Z_i,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;}]\&Proportional;\mathrm{Gamma}\{a+\frac{1}{2},b+\frac{1}{2}\mathrm{\&gamma;}{\left(K\left(X\right)\right)}^{-1}\mathrm{\&gamma;}\}$

由于所有参数的满条件分布均可得到闭合形式，因此使用Gibbs抽样算法循环计算各个参数的后验分布。Gibbs抽样算法如下：

步骤1：初始化参数[Z⁽⁰⁾，α⁽⁰⁾，γ⁽⁰⁾，τ⁽⁰⁾]。

步骤2：在第t次迭代时，

(i)求Z^(t)：[Z|Z^(t-1)，α^(t-1)，γ^(t-1)]；

(ii)求α^(t)：[α|Z^(t)，α^(t-1)，γ^(t-1)]；

(iii)求γ^(t)：[γ|Z^(t)，α^(t)，τ^(t-1)]；

(iv)求τ^(t)：[τ|γ^(t)]。

步骤3：增加迭代次数t直至到达规定的迭代次数。

通过上述Gibbs抽样算法，可得到马尔可夫链蒙特卡洛(MCMC)样本用于贝叶斯推断和预测，计算95％的贝叶斯置信区间，并得到点估计量，由此得到τ^-1，检验τ^-1是否为0或τ^-1K(X)是否等于τ^-1I。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (6)

1.一种图像分析方法，其特征在于，包括如下步骤：

获取一组T1加权图像，包括对照组图像和实验组图像；

对所有所述T1加权图像预处理；

将预处理后的T1加权图像非线性配准到预设的模板图像中；

将非线性配准后的T1加权图像进行分割，产生多个感兴趣区域；

提取每个所述感兴趣区域的至少一种特征数据，得到特征矩阵信息；

根据预设的贝叶斯多层模型和所述特征矩阵信息，计算得到一参数值，所述预设的贝叶斯多层模型如下：

Pr{Y＝1|C，γ(X)}＝Φ{α^TC+γ(X)}，

γ(X)～MN{0，τ^-1K(X)}，

其中，Y_i为二进制响应变量，Y_i＝1表示第i幅图像来自实验组，Y_i＝0表示第i幅图像来自对照组，i∈{1，2，…，n}，Y_n×1为二进制响应结果，C_n×1为协变量，X_n×p为所述特征矩阵信息，Φ(.)为标准正态累积分布函数，α为回归系数，γ(X)＝{γ(x₁)，γ(x₂)，…，γ(x_n)}表示均值为0，协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)的高斯随机过程，τ服从Gamma分布，K是一个n×n矩阵，K_ij＝K(x_i，x_j)，K为核函数，所述参数值为所述预设的贝叶斯多层模型中τ的值；

根据所述参数值与预设阀值的关系，判断所述实验组图像和所述对照组图像之间感兴趣区域的差异显著性。

2.根据权利要求1所述的，其特征在于，所述提取每个感兴趣区域的至少一种特征数据，得到特征矩阵信息步骤是从所述感兴趣区域R_q，q∈{1，2，…，Q}中分别提取每幅图像的多种特征信息：

其中i为图像索引号i∈{1，2，…，n}，q为感兴趣区域索引号q∈{1，2，…，Q}，p为该区域特征信息的维度p≥1，对每个感兴趣区域R_q得到的所述特征矩阵信息如下：

其中x_ij表示第i幅图像感兴趣区域的第j个参数元素，i∈{1，2，…，n}，j∈{1，2，…，p}。

3.根据权利要求1所述的，其特征在于，所述预设的贝叶斯多层模型计算方式如下：

首先，定义了一个隐式变量Z＝(z₁，z₂，…，z_n)，Z～MN{α^TC+γ(X)，I}，使得

$Y_i=\left\{\begin{array}{cc}1& \mathrm{if}{Z}_i\&GreaterEqual;0\\ 0& \mathrm{if}{Z}_i<0\end{array}\right.,$

将协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)中的K(x，x′)定义为高斯核，则：

$K(x,x^{\&prime;})=e^{-\frac{\left|\right|x-x^{\&prime;}\left|\right|}{2}},$

然后，假设γ(X)～MN(0，τ^-1K(X)}，α～MN(0，φ^-1)和τ～Gamma(a，b)，则联合后验分布为：

$[Z,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;},\mathrm{\&tau;}|y,C,X]\&Proportional;\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}(I_{Z_i\&GreaterEqual;0}{I}_{y_i=1}+I_{Z_i<0}{I}_{y_i=0})\&times;\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1)$

$\&times;\mathrm{MN}\{\mathrm{\&gamma;};0,{\mathrm{\&tau;}}^{-1}K\left(X\right)\}\&times;\mathrm{Gamma}(\mathrm{\&tau;};a,b)$

$\&times;N({\mathrm{\&alpha;}}_0;0,{\mathrm{\&phi;}}^{-1})\&times;N({\mathrm{\&alpha;}}_1;0,{\mathrm{\&phi;}}^{-1}),$

得到Z_i的满条件概率分布的闭合形式：

$\left[Z_i\right|y_i=1,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,{\mathrm{\&gamma;}}_i,\mathrm{\&tau;}]\&Proportional;I_{Z_i<0}{I}_{y_i=0}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1)$

$\left[Z_i\right|y_i=0,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,{\mathrm{\&gamma;}}_i,\mathrm{\&tau;}]\&Proportional;I_{Z_i\&GreaterEqual;0}{I}_{y_i=1}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1),$

该分布为截断正态分布。α，γ，τ的满条件分布求得如下：

$\left[{\mathrm{\&alpha;}}_0\right|Z_i,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;},\mathrm{\&tau;}]\&Proportional;N\{\frac{{\mathrm{\&Sigma;}}_{i=1}^n(Z_i-{\mathrm{\&alpha;}}_1{C}_i-{\mathrm{\&gamma;}}_i)}{n+\mathrm{\&phi;}},{(n+\mathrm{\&phi;})}^{-1}\}$

$\left[{\mathrm{\&alpha;}}_1\right|Z_i,{\mathrm{\&alpha;}}_0,\mathrm{\&gamma;},\mathrm{\&tau;}]\&Proportional;N\{\frac{{\mathrm{\&Sigma;}}_{i=1}^n{C}_i(Z_i-{\mathrm{\&alpha;}}_0-{\mathrm{\&gamma;}}_i)}{{\mathrm{\&Sigma;}}_{i=1}^n{C}_i^2+\mathrm{\&phi;}},{(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{C}_i^2+\mathrm{\&phi;})}^{-1}\}$

[γ|Z_i，α₀，α₁，τ]∝MN{[1+(τ^-1K(X))^-1](Z-α′C)，[1+τ^-1K(X))^-1]^-1}

$\left[\mathrm{\&tau;}\right|Z_i,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;}]\&Proportional;\mathrm{Gamma}\{a+\frac{1}{2},b+\frac{1}{2}\mathrm{\&gamma;}{\left(K\left(X\right)\right)}^{-1}\mathrm{\&gamma;}\};$

使用Gibbs抽样算法循环计算各个参数的后验分布，Gibbs抽样算法如下：

步骤1：初始化参数[Z⁽⁰⁾，α⁽⁰⁾，γ⁽⁰⁾，τ⁽⁰⁾]；

步骤2：在第t次迭代时，

(i)求Z^(t)：[Z|Z^(t-1)，α^(t-1)，γ^(t-1)]；

(ii)求α^(t)：[α|Z^(t)，α^(t-1)，γ^(t-1)]；

(iii)求γ^(t)：[γ|Z^(t)，α^(t)，τ^(t-1)]；

(iv)求τ^(t)：[τ|γ^(t)]。

步骤3：增加迭代次数t直至到达规定的迭代次数；

通过上述Gibbs抽样算法，得到马尔可夫链蒙特卡洛(MCMC)样本用于贝叶斯推断和预测，计算95％的贝叶斯置信区间并得到点估计量，由此得到τ^-1。

4.一种图像分析系统，其特征在于，包括获取单元、处理单元、配准单元、分割单元、提取单元、计算单元和分析单元，

所述获取单元用于获取一组T1加权图像，包括对照组图像和实验组图像；

所述处理单元用于对所有所述T1加权图像预处理；

所述配准单元用于将预处理后的T1加权图像非线性配准到预设的模板图像中；

所述分割单元用于将非线性配准后的T1加权图像进行分割，产生多个感兴趣区域；

所述提取单元用于提取每个所述感兴趣区域的至少一种特征数据，得到特征矩阵信息；

所述计算单元用于根据预设的贝叶斯多层模型和所述特征矩阵信息，计算得到一参数值，所述预设的贝叶斯多层模型如下：

Pr{Y＝1|C，γ(X)}＝Φ{α^TC+γ(X)}，

γ(X)～MN{0，τ^-1K(X)}，

其中，Y_i为二进制响应变量，Y_i＝1表示第i幅图像来自实验组，Y_i＝0表示第i幅图像来自对照组，i∈{1，2，…，n}，Y_n×1为二进制响应结果，C_n×1为协变量，X_n×p为所述特征矩阵信息，Φ(.)为标准正态累积分布函数，α为回归系数，γ(X)＝{γ(x₁)，γ(x₂)，…，γ(x_n)}表示均值为0，协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)的高斯随机过程，τ服从Gamma分布，K是一个n×n矩阵，K_ij＝K(x_i，x_j)，K为核函数，所述参数值为所述预设的贝叶斯多层模型中τ的值；

所述分析单元用于根据所述参数值与预设阀值的关系，判断所述实验组图像和所述对照组图像之间感兴趣区域的差异显著性。

5.根据权利要求4所述的，其特征在于，所述提取每个感兴趣区域的至少一种特征数据，得到特征矩阵信息步骤是从感兴趣区域R_q，q∈{1，2，…，Q}中分别提取每幅图像的多种特征信息：

其中i为图像索引号i∈{1，2，…，n}，q为感兴趣区域索引号q∈{1，2，…，Q}，p为该区域特征信息的维度p≥1，对每个感兴趣区域R_q得到的所述特征矩阵信息如下：

其中x_ij表示第i幅图像感兴趣区域的第j个参数元素，i∈{1，2，…，n}，j∈{1，2，…，p}。

6.根据权利要求4所述的，其特征在于，所述预设的贝叶斯多层模型计算方式如下：

首先，定义了一个隐式变量Z＝(z₁，z₂，…，z_n)，Z～MN{α^TC+γ(X)，I}，使得

$Y_i=\left\{\begin{array}{cc}1& \mathrm{if}{Z}_i\&GreaterEqual;0\\ 0& \mathrm{if}{Z}_i<0\end{array}\right.,$

将协方差cov{γ(x)，γ(x′)}＝τ^-1K(x，x′)中的K(x，x′)定义为高斯核，则：

$K(x,x^{\&prime;})=e^{-\frac{\left|\right|x-x^{\&prime;}\left|\right|}{2}},$

然后，假设γ(X)～MN{0，τ^-1K(X)}，α～MN(0，φ^-1)和τ～Gamma(a，b)，则联合后验分布为：

$[Z,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;},\mathrm{\&tau;}|y,C,X]\&Proportional;\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}(I_{Z_i\&GreaterEqual;0}{I}_{y_i=1}+I_{Z_i<0}{I}_{y_i=0})\&times;\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1)$

$\&times;\mathrm{MN}\{\mathrm{\&gamma;};0,{\mathrm{\&tau;}}^{-1}K\left(X\right)\}\&times;\mathrm{Gamma}(\mathrm{\&tau;};a,b)$

$\&times;N({\mathrm{\&alpha;}}_0;0,{\mathrm{\&phi;}}^{-1})\&times;N({\mathrm{\&alpha;}}_1;0,{\mathrm{\&phi;}}^{-1}),$

得到Z_i的满条件概率分布的闭合形式：

$\left[Z_i\right|y_i=1,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,{\mathrm{\&gamma;}}_i,\mathrm{\&tau;}]\&Proportional;I_{Z_i<0}{I}_{y_i=0}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1)$

$\left[Z_i\right|y_i=0,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,{\mathrm{\&gamma;}}_i,\mathrm{\&tau;}]\&Proportional;I_{Z_i\&GreaterEqual;0}{I}_{y_i=1}N(Z_i;{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_1{C}_i+{\mathrm{\&gamma;}}_i,1),$

该分布为截断正态分布。α，γ，τ的满条件分布求得如下：

$\left[{\mathrm{\&alpha;}}_0\right|Z_i,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;},\mathrm{\&tau;}]\&Proportional;N\{\frac{{\mathrm{\&Sigma;}}_{i=1}^n(Z_i-{\mathrm{\&alpha;}}_1{C}_i-{\mathrm{\&gamma;}}_i)}{n+\mathrm{\&phi;}},{(n+\mathrm{\&phi;})}^{-1}\}$

$\left[{\mathrm{\&alpha;}}_1\right|Z_i,{\mathrm{\&alpha;}}_0,\mathrm{\&gamma;},\mathrm{\&tau;}]\&Proportional;N\{\frac{{\mathrm{\&Sigma;}}_{i=1}^n{C}_i(Z_i-{\mathrm{\&alpha;}}_0-{\mathrm{\&gamma;}}_i)}{{\mathrm{\&Sigma;}}_{i=1}^n{C}_i^2+\mathrm{\&phi;}},{(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{C}_i^2+\mathrm{\&phi;})}^{-1}\}$

[γ|Z_i，α₀，α₁，τ]∝MN{[1+(τ^-1K(X))^-1](Z-α′C)，[1+τ^-1K(X))^-1]^-1}

$\left[\mathrm{\&tau;}\right|Z_i,{\mathrm{\&alpha;}}_0,{\mathrm{\&alpha;}}_1,\mathrm{\&gamma;}]\&Proportional;\mathrm{Gamma}\{a+\frac{1}{2},b+\frac{1}{2}\mathrm{\&gamma;}{\left(K\left(X\right)\right)}^{-1}\mathrm{\&gamma;}\};$

使用Gibbs抽样算法循环计算各个参数的后验分布，Gibbs抽样算法如下：

步骤1：初始化参数[Z⁽⁰⁾，α⁽⁰⁾，γ⁽⁰⁾，τ⁽⁰⁾]；

步骤2：在第t次迭代时，

(i)求Z^(t)：[Z|^(t-1)，α^(t-1)，γ^(t-1)]；

(ii)求α^(t)：[α|Z^(t)，α^(t-1)，γ^(t-1)]；

(iii)求γ^(t)：[γ|Z^(t)，α^(t)，τ^(t-1)]；

(iv)求τ^(t)：[τ|γ^(t)]。

步骤3：增加迭代次数t直至到达规定的迭代次数；

通过上述Gibbs抽样算法，得到马尔可夫链蒙特卡洛(MCMC)样本用于贝叶斯推断和预测，计算95％的贝叶斯置信区间并得到点估计量，由此得到τ^-1。

Images