对数正态多层贝叶斯模型的参数估计.pdf

第 35 卷第 3 期湖南文理学院学报(自然科学版)Vol.35 No.32023 年 9 月Journal of Hunan University of Arts and Science(Science and Technology)Sep.2023doi:10.3969/j.issn.16726146.2023.03.003对数正态多层贝叶斯模型的参数估计 王志凯,黄介武(贵州民族大学 数据科学与信息工程学院,贵州 贵阳,550025)摘要:为了更好地捕捉呈偏态分布数据的变化,提高统计推断的精确度,将对数正态多层先验分布的构造方法与贝叶斯定理结合建立了对数正态多层贝叶斯模型。利用Gibbs抽样算法对各未知参数进行贝叶斯估计,并对使用 Gibbs 算法所生成的迭代链进行收敛性诊断。随机模拟结果显示,在相对误差、均方误差(MSE)准则下,贝叶斯估计的效果较似然估计更优。最后,通过实证分析证明了所建立的模型是切实可行的。关键词:正态多层贝叶斯模型;对数正态多层贝叶斯模型;贝叶斯估计;Gibbs 算法中图分类号:O 212.8文献标志码:A文章编号:16726146(2023)03001208Parameter estimation of lognormal multilayer Bayesian modelWang Zhikai,Huang Jiewu(School of Data Science and Information Engineering,Guizhou Minzu University,Guiyang 550025,China)Abstract:In order to better capture the changes of skewed distribution data and improve the accuracy of statisticalinference,this paper combines the construction method of lognormal multilayer prior distribution with Bayesiantheorem to establish a lognormal multilayer Bayesian model.The Gibbs sampling algorithm is used to estimate theunknown parameters,and the convergence of the iterative chain generated by the Gibbs algorithm is diagnosed.Therandom simulation results show that the Bayesian estimation is better than the maximum likelihood estimation underthe relative error and MSE(mean square error)criteria.Finally,the empirical analysis proves that the establishedmodel is feasible.Key words:normal multilayer Bayesian model;lognormal multilayer Bayesian model;Bayesian estimation;Gibbsalgorithm多层贝叶斯模型是具有结构化层次的统计模型,它可以用来为复杂的统计问题建立多层模型从而避免参数过多导致的过拟合问题,其应用非常广泛,特别是在水文统计12、数据处理34、金融经济预测、地震强度预测等领域。在现实生活中,服从或近似于正态分布的数据具有不同的层次结构。针对不存在组内和组间异质性的数据,利用正态单层贝叶斯模型对其进行统计推断得到的结果较好。针对具有嵌套和分层结构的数据,若再用正态单层贝叶斯模型来拟合这些数据,会出现过多的参数导致过拟合的问题。因此,AndrewGelman5为了解决这些数据过拟合的问题,首次提出两层正态贝叶斯模型,并讨论了均值参数、方差参数的先验分布选择问题。此外,国内外学者还进一步对正态多层贝叶斯模型进行了很多研究。例如,刘金山等67将 Andrew Gelman 提出的两层正态贝叶斯模型进行扩充到多层正态贝叶斯模型,并基于正态多层混合贝叶斯模型研究大脑 FMRI 图像分割;针对小样本情况,Yang 等8用极大似然方法和贝叶斯方通信作者:黄介武,。收稿日期:20221026基金项目:贵州省科技计划基金项目(黔科合基础20171083 号);贵州省基础研究计划(软科学)(黔科合支201920001)。第 3 期王志凯,等:对数正态多层贝叶斯模型的参数估计13法对正态多层贝叶斯模型中的未知参数进行了估计,并用迭代轨迹图对所得的贝叶斯估计进行了收敛性诊断;朱彬彬9基于 MCMC 算法对正态多层贝叶斯模型中的未知参数进行了估计,使用图像判别、遍历均值和方差法对所得的贝叶斯估计进行了收敛性诊断;Haydar Demirhan 等10研究了正态多层贝叶斯模型中方差参数的联合先验分布。上面提到的研究都是在数据服从正态分布的假设下进行讨论的,而实际的数据呈现严格对称分布的情况非常的少,特别是经济学、气候科学和生态学等领域中的实际数据常常会呈现偏态分布。此时,如果再继续基于正态假设下用这些数据作统计推断,可能会获得不合理甚至错误的结论。为了更好的捕捉到这些呈偏态分布数据的变化,进而提高统计推断的精确度。常用的数据处理方法是对这些数据进行取对数,这种数据处理方法可以在一定程度上修正数据的偏态分布问题,使其更加接近于正态分布。因此,本文为了解决实际中一些数据呈偏态分布而不适用于正态多层贝叶斯模型的问题,将对数正态多层先验分布的构造方法与贝叶斯定理结合,建立对数正态多层贝叶斯模型,利用贝叶斯方法对对数正态多层贝叶斯模型进行参数估计,并通过实证分析说明该模型能更好的对呈偏态分布的数据进行统计推断。1 模型1.1多层贝叶斯模型多层贝叶斯方法是一种以贝叶斯原理为基础,应用多层先验分布法确定各层次的先验分布,再进行统计推断的贝叶斯建模方法。多层贝叶斯方法最核心的思想就是将多层先验分布的构造方法与贝叶斯定理进行结合,从而完成对多层贝叶斯模型的构建。假设=(1,2,k)为一组待估计的未知参数,是的参数空间,X=(X1,X2,Xk)为通过观察或试验所获得的分层样本数据集,其中12(,),1,2,iiiiinXxxxik。如果在观察这些分层样本数据之前,根据以往的经验或历史资料获得参数的总体分布信息,且该总体分布需要通过引入与相关的另外一个未知参数来描述。那么对多层贝叶斯模型的构建包括以下几个步骤:(1)构建的先验分布12(,|)k ,其中为超参数;(2)利用其他可利用的信息 h 构建超参数的先验分布(|)h,其中 h 为已知信息;(3)通过贝叶斯定理,将先验分布与似然函数相乘,求得模型参数的后验分布。通过以上描述的构造过程,即得到多层贝叶斯模型:总体分布|(),xp x;xX第一层先验|1(),参数;第二层先验2(|),h超参数,h为已知信息。1.2正态多层贝叶斯模型给定一个分层数据集 X=(X1,X2,Xn),其中12(,),(1,2,)iiiiinXxxxim,假设该分层数据集满足多层贝叶斯模型:2,221122022233(,),1,2,1,2,;(,),(,),1,2,;(,),(,);(,)i jiiixNjn imNIGa a bimNIGa a bIGa a b (1)则称该模型为正态多层贝叶斯模型。其中,(,)iiIGa a b表示参数为(,)iia b的逆伽马(IGamma)分布,0,ia(1,2,3)ib i 是给定的超参数。1.3对数正态多层贝叶斯模型为了解决实际中一些数据呈偏态分布而不适用于正态多层贝叶斯模型的问题,本文将对数正态多层先验分布的构造方法与贝叶斯定理结合,建立如下所示的对数正态多层贝叶斯模型。给定一个分层数据集 Y=(Y1,Y2,Yk),其中12(ln,ln,ln),iiiiinYxxx假设该分层数据集满足以下的多层贝叶斯模型:14湖南文理学院学报(自然科学版)2023 年222112200222033ln(,),1,2,1,2,;(,),(,),1,2,;(,),(,);(,)ijiiixNjn ikNIGaikNIGaIGa (2)则称该模型为对数正态多层贝叶斯模型。其中,(,)iiIGa 表示参数为(,)ii 的逆伽马(IGamma)分布,2220,i 是未知的参数,0,1,2,3iii 是给定的超参数。2 模型的参数估计2.1对数正态多层贝叶斯模型参数的极大似然估计设12(,)kYY YY是来自模型(2)的观测数据,则有似然函数为2222/22/2200222221100()11111(|,)()exp(ln)()exp()exp()222inknkijiijiL Yx 31211131222222200111()exp()()exp()()exp(),1,2,ik(3)其中:12(ln,ln,ln)iiiiinYxxx,12(,)k,1kiinn。对式(3)两边取对数,可得样本的对数似然函数为222220222221100111111ln(|,)ln()(ln)ln()()ln()(22222inkijiijinkL Yx 2312012322222200111)(1)ln()(1)ln()(1)ln(),1,2,ik。(4)再根据式(4),在固定时,对i,2,2,20求极大似然估计,得2201111,11(ln)/(1),()/(1),(ln)/(1)22inkiijiiijijii jnxnkx,(5)22222200331111()/(1),()/(1),2222kiik(6)其中,1,2,ik。2.2对数正态多层贝叶斯模型参数的贝叶斯估计由模型(2)可知,设12(,)kYY YY是来自对数正态分布2(,)N 的观测数据,其中12(,12),(ln,ln,ln)ikiiiinYxxx,1,2,ik。此样本的联合密度函数为2222/2220,11(|,)()exp(ln)/2),1,2,knijiii jip Yxik nn 。(7)得到参数2220,的条件后验密度分别为:(1)22222200(|,)(|,)(|iiYp Y 22222222222222222111,)exp(ln)exp()exp(),222iniiijiiiijiinnxYnn 1,2,ik其中,111,ln,inkiiijijinnYxn2(|,)i 是i的先验分布。(2)2220(|,)(Y 2222222200000222222221000011|,)(|)exp()exp()exp(222kiikkkk 2),其中,1()/,kiik20(|)是的先验分布。(3)22222200(|,)(|,)(Yp Y 第 3 期王志凯,等:对数正态多层贝叶斯模型的参数估计15111/212/22221122222,1,11111)()exp(ln)()exp()()exp(ln)/(2),2nnijiijii ji jxx其中,2()是2的先验分布。(4)22222202111(|,)(|,)()exp()2kiiY 221/212222222111()exp()()exp()/(2),kkii其中,2()是2的先验分布。(5)20(|33122112222030320022222000000()()111,)(|)()exp()()exp()()exp(),22Y 其中,20()是20的先验分布。进而得到各参数的条件后验分布为:222222220222222(|,)(,),iiiiiinYNYnnn 1,2,;ik22222220000222222000(|,)(,)kYNkkk ;2220(|,)(2nYIGa 211,11,(ln)2ijii jx;222202211(|,)(,()22kiikYIGa ;2220(|,)Y 2303(1/2,()/2)IGa。由于所得到的各未知参数的条件后验分布为标准分布,可直接用Gibbs抽样进行抽样,得到每个未知参数条件后验分布的平稳分布。Gibbs 抽样是 MCMC 方法里最简单且应用最广泛的抽样方法11,则Gibbs 抽样算法的具体步骤如下:(1)给定目标平稳分布2220(,)(1,2,)iik ,假设需要进行N次迭代得到 Gibbs样本量为N,在第 A 次迭代时的 Gibbs 抽样开始收敛;(2)给定初始值为(0)(0)2(0)2(0)2(0)0,1,2,iik;(3)对于1,2,tAN,假设第1t 次迭代的估计值为(1)(1)2(1)2(1)2(1)0(,)ttttti,则第t次迭代的步骤为:(a)从条件后验分布(1)2(1)2(1)2(1)0(|,)ttttiY (1,2,)ik随机抽取()ti,用()ti更新i;(b)从条件后验分布(1)2(1)2(1)2(1)0(|,)ttttY 随机抽取()t,用()t更新;(c)从条件后验分布2(1)(|,ti(1)2(1)2(1)0,)tttY随机抽取2()t,用2()t更新2;(d)从条件后验分布2(1)(1)2(1)2(1)0(|,)ttttiY 随机抽取2()t,用2()t更新2;(e)从条件后验分布2(1)(1)0(|,tti 2(1)2(1),)ttY随机抽取2()0t,用2()0t更新20;(4)重复以上第(3)过程,迭代 N 次,得到(1)(1)2(1)2(1)2(1)(2)(2)2(2)2(2)2(2)00(,),(,),ii()()2()2()2()()()2()2()2()00(,),(,),AAAAANNNNNii1,2,ik,即为根据 Gibbs 抽样获得的样本集。假设得到各参数的Gibbs样本容量为N,且在第A次迭代时开始收敛,则把后面NA个样本的均值作为各未知参数的贝叶斯估计值。即()()22()2111111,NNNtttiitAtAtANANANA2()22()001111,NNtttAtANANA其中,1,2,ik。3 模拟研究基于上述对对数正态多层贝叶斯模型参数的极大似然估计和贝叶斯估计的理论研究,下面利用 R软件进行随机模拟。分别取样本量 n1=n2=50,n1=n2=100,n1=n2=200,0=5,i=i=3(i=1,2,3),真实值(1,2,2,2,20)=(3.4,3.8,4.1,0.4,1.2,1.5)时,根据极大似然估计式(5)、(6),对未知参数1,2,2,2,16湖南文理学院学报(自然科学版)2023 年20进行极大似然估计,再根据各未知参数的条件后验分布利用 Gibbs 算法进行随机模拟,把所得到的样本均值作为未知参数的最终贝叶斯估计值。在模拟的过程中,先进行 1 104次的预迭代以保证收敛,在预迭代的基础上再进行 2 104次的迭代,最后将 2 104个样本的均值作为各未知参数的贝叶斯估计。两种估计方法在不同样本量中的模拟结果分别如表 13 所示。表 13 给出了在不同样本量下,用两种估计方法求得的各未知参数估计值的情况。只有当运用Gibbs 算法最终生成的迭代链收敛时,才能证明求得的模型中各未知参数的贝叶斯估计值是可靠的。因此,利用图像法对 Gibbs 抽样生成的迭代链进行收敛性诊断。在样本量不同的情况下,对该模型中各未知参数的 Gibbs 抽样所生成的迭代链进行模拟轨迹图,如图 13 所示。表 1n1=n2=50 各未知参数的估计值、相对误差与 MSE参数真实值估计值相对误差MSEMLE贝叶斯估计MLE贝叶斯估计MLE贝叶斯估计13.43.280 33.432 20.035 10.009 42.86 e041.03e0523.83.557 33.712 40.063 80.023 01.17e047.66e054.13.945 94.031 40.037 50.016 72.37 e044.69e0520.40.367 80.382 50.080 40.023 61.03e053.04e0621.21.087 51.224 90.093 70.020 71.02 e046.22e06201.50.989 81.468 50.334 10.020 90.002 39.86e06表 2n1=n2=100 各未知参数的估计值、相对误差与 MSE参数真实值估计值相对误差MSEMLE贝叶斯估计MLE贝叶斯估计MLE贝叶斯估计13.43.442 53.421 30.016 30.006 23.10e052.28e0623.83.668 53.745 90.034 50.014 21.72 e041.45e054.14.004 24.031 00.023 30.016 84.58e052.37e0520.40.391 10.392 90.022 10.017 53.92e072.46e0721.21.091 31.224 50.090 50.020 45.89e051.50e06201.50.998 71.469 20.340 00.020 41.25 e042.36e06表 3n1=n2=200 各未知参数的估计值、相对误差与 MSE参数真实值估计值相对误差MSEMLE贝叶斯估计MLE贝叶斯估计MLE贝叶斯估计13.43.379 83.418 70.005 90.005 52.02e068.75e0723.83.728 33.768 20.018 80.008 32.56e052.52e064.14.036 04.049 30.015 50.012 31.02e056.43e0620.40.395 00.396 00.012 20.009 96.03e083.94e0821.21.098 11.215 70.084 10.013 13.16e051.24e06201.51.015 81.479 10.322 70.013 90.65 e042.18e06051203.63.21510迭代次数/103迭代次数/103迭代次数/1034.023.73.405201510052015101357图 1各未知参数的 Gibbs 抽样迭代轨迹图(n1=n2=50)051020150.50.3迭代次数/103迭代次数/103迭代次数/103051020150510201501020001020第 3 期王志凯,等:对数正态多层贝叶斯模型的参数估计17由图 13 可知,所有的 Gibbs 抽样迭代轨迹图没有明显的周期和趋势,说明在样本量不同的情况下,模型中各未知参数通过 Gibbs 抽样生成的所有迭代链都是收敛的。因此,将后面 2 104个样本的均值作为各未知参数的贝叶斯估计值是合理的。表 13 的模拟结果表明,样本量 n 相同时,各未知参数的贝叶斯估计值都比极大似然估计值更接近于真实值。在相对误差准则下,贝叶斯估计的相对误差均不超过 3%,极大似然估计的相对误差均不超过 9%,说明了贝叶斯估计的精确度高于极大似然估计;在 MSE 准则下,贝叶斯估计的效果优于极大似然估计,说明了 Gibbs 抽样得到的迭代值的波动较小,相对稳定。针对样本量 n 不同的情况,当样本量 n 逐渐增大时,参数估计的相对误差、MSE 逐渐减小,说明了随着样本量 n 的增加,贝叶斯估计的精度越来越高,Gibbs抽样得到的迭代值的波动越小。4 实证分析本文实证分析的数据是选自文献12记录对 20 个实验材料进行测试得到的相应最大应力数据,测试记录的数据如表 4 所示。4.1数据的正态性检验首先,对选取的两组数据分别取对数,再采用Shapiro-Wilk检验对选取的数据进行正态性检验。假设检验分别为:H0(这两组数据服从对数正态分布);H1(这两组数据不服从对数正态分布)。其中显著性水平取值为 0.05。利用 Shapiro-Wilk 检验方法进行检验得到的结果如表 5 所示。由表 5 可知,这两组数据的 Shapiro-Wilk检验的 P 值都是大于 0.05，说明原假设成立,即这两组数据服从对数正态分布。因此,可以用这些数据来对本文所提出的模型中的未知参数进行估计。4.2模型参数的估计在模型(2)中,x1,x2分别是记录两组实验材料测试的最大应力,k=2,n1=n2=10,给定超参数0=6,i=i=2(i=1,2,3)。根据极大似然估计的式(5)、(6),对未知参数1,2,2,2,图 2各未知参数的 Gibbs 抽样迭代轨迹图(n1=n2=100)051203.53.21510迭代次数/103迭代次数/103迭代次数/1034.023.83.605201510052015101357051020150.500.30迭代次数/103迭代次数/103迭代次数/1030510201505102015010200020图 3各未知参数的 Gibbs 抽样迭代轨迹图(n1=n2=200)051203.453.251510迭代次数/103迭代次数/103迭代次数/10323.83.60520151005201510246051020150.450.30迭代次数/103迭代次数/103迭代次数/1030510201505102015010150030515表 4测试的数据x1x2448.1458.7449.4456.1450.8448.1445.6452.1457.4450.8452.1442.8448.1449.4444.1445.5450.8450.8456.1452.1表 5Shapiro-Wilk检测结果数据组P 值x10.767 6x20.945 518湖南文理学院学报(自然科学版)2023 年20进行极大似然估计,再根据各未知参数的条件后验分布,利用 Gibbs 抽样算法来估计模型中的未知参数。假设在抽样的过程中,先进行 1 104次的预迭代以保证收敛,在预迭代的基础上再进行 2 104次的迭代,最终得到模型中各未知参数的估计值如表 6 所示。表 6 给出了进行 3 104次迭代时,用两种估计方法求得模型的各未知参数估计值的情况。只有当 Gibbs 算法最终生成的迭代链收敛时,才能证明求得的模型中各未知参数的贝叶斯估计值是可靠的。因此,利用方差比法对 Gibbs 抽样生成的迭代链进行收敛性诊断。结合上述的实际数据以及构建的对数正态多层贝叶斯模型,选取 6 个未知参数在 Gibbs 抽样下所生成的 6 条迭代链,每条迭代链的迭代次数为 1 104次,最终可以计算得出模型中各个未知参数的潜尺度减缩因子(potential scale reductionfactor)值。若潜尺度减缩因子的值1.1,则可以认为此时 Gibbs 抽样所生成的迭代链是收敛的。各未知参数潜尺度减缩因子的估计值如表 7 所示(表 7 中的 Point est.表示模型中各未知参数潜尺度减缩因子的点估计,Upper C.I.表示模型中各未知参数潜尺度减缩因子估计值的置信水平为 97.5%的置信区间的上限)。只有当模型中某个未知参数潜尺度减缩因子估计值的置信区间的上限接近于 1,才能说明该参数在 Gibbs 抽样下所生成的迭代链是收敛的。即证明了求得该参数的贝叶斯估计值是可靠的。由表 7 可知,模型中各未知参数潜尺度减缩因子估计值的置信区间的上限都接近于1,证明了该模型中各未知参数在 Gibbs 抽样下生成的所有迭代链都是收敛的。即该模型中各未知参数的贝叶斯估计值是可靠的。为了避免模型中某个未知参数潜尺度减缩因子的 97.5%置信区间的上限偶尔接近于 1,而导致这个参数的 Gibbs 抽样迭代链被误判定为已经收敛的问题,给出了该模型中各未知参数的潜尺度减缩因子的演化图,如图 4 所示。通过图 4 可以更加准确的判断各未知参数的潜尺度减缩因子是否真的收敛或者还存在震荡变化。表 6模型各未知参数的估计值估计方法估计值122220MLE5.181 6005.182 385.454 660.965 121.037 170.859 478贝叶斯估计6.109 3256.109 366.078 840.199 721.167 841.442 430表 7收敛性诊断结果参数Point est.Upper C.I.11.001.0021.011.011.001.0021.001.0021.001.00201.011.01迭代次数/10302468101.01.21.41.6潜尺度减缩因子1median97.5%迭代次数/10302468101.001.101.20潜尺度减缩因子2median97.5%迭代次数/10302468101.001.10潜尺度减缩因子median97.5%迭代次数/10302468101.01.52.02.5潜尺度减缩因子median97.5%第 3 期王志凯,等:对数正态多层贝叶斯模型的参数估计19由图 4 可知,各未知参数均在 1 104这个时点之后的迭代链才是真正的收敛。说明了要先进行 1 104次的预迭代才能保证得到的各个未知参数贝叶斯估计值是可靠的。5 总结本文为了解决实际中一些数据呈偏态分布在正态多层贝叶斯模型下进行统计推断得到的结果不理想的问题,进而提出了对数正态多层贝叶斯模型,并对所提出模型的未知参数进行了极大似然估计和贝叶斯估计。随机模拟结果表明,相较于极大似然估计,各未知参数的贝叶斯估计值更稳定,也更接近于真实值。并通过实证分析证明了所提出的模型是切实可行的。参考文献:1BRACKEN C,RAJAGOPALANR B,CHENG L,et al.Spatial Bayesian hierarchical modeling of precipitation extremesover a large domain J.Water Resources Research,2016,52(8):6 6436 655.2YAN H X,HAMID M.A regional Bayesian hierarchical model for flood frequency analysis J.Stochastic EnvironmentalResearch and RiskAssessment,2015,29(3):1 0191 036.3文雯,文小焱,胡珊,等.贝叶斯层次模型在嵌套结构调查数据中的应用研究J.中国卫生统计,2015,32(2):190193.4张天嵩.基于二项式-正态层次模型框架下比例的贝叶斯 Meta 分析方法及实现J.中国循证儿科杂志,2019,14(2):123128.5GELMAN A.Prior distributions for variance parameters in hierarchical models(comment on article by Browne andDraper)J.Bayesian analysis,2006,1(3):515534.6刘金山,夏强.基于 MCMC 算法的贝叶斯统计方法M.北京:科学出版社,2016.7刘金山,陈镇坤,黄来华.基于贝叶斯分层混合模型的大脑 FMRI 图像分割J.数理统计与管理,2015,34(4):603611.8YANG J W,WANG J H,ZHAO C Q.Study on New Hierarchical Bayesian Model Based on Normal Distribution for SmallSample Size J.Journal of Computational and Theoretical Nanoscience,2016,13(2):1 3231 329.9朱彬彬.基于 MCMC 的贝叶斯统计方法在分层正态模型中的应用D.兰州:兰州大学,2018.10 HAYDAR D,ZEYNEP K.Joint prior distributions for variance parameters in Bayesian analysis of normal hierarchicalmodels J.Journal of MultivariateAnalysis,2015,135:163174.11 韦学永,金良琼,沈婷,等.一种基于 ALD 的贝叶斯分位数回归估计J.湖南文理学院学报(自然科学版),2022,34(2):17.12 POBOKOV Ivana,SEDLIAKOV Zuzana,MICHALKOV Mria,et al.Monte Carlo comparison of estimationmethods for the two-parameter lognormal distribution J/OL.MATEC Web of Conferences,2019,254(205):08002.https:/doi.org/10.1051/matecconf/201925408002.(责任编校:张红)迭代次数/10302468101.01.41.82.2潜尺度减缩因子median97.5%迭代次数/10302468101.01.41.82.2潜尺度减缩因子0median97.5%图 4各参数的潜尺度减缩因子演化图