拉普拉斯分布参数的近似贝叶斯估计.pdf

应用概率统计第 40 卷第 1 期2024 年 2 月Chinese Journal of Applied Probability and StatisticsFeb.,2024,Vol.40,No.1,pp.18-32doi:10.3969/j.issn.1001-4268.2024.01.002拉普拉斯分布参数的近似贝叶斯估计杨彦娇王立春(北京交通大学数学与统计学院,北京,100044)摘要:拉普拉斯分布是刻画尖峰厚尾数据的重要分布之一.本文提出拉普拉斯分布两参数具有显式解的线性近似贝叶斯估计,通过理论证明和数值模拟验证了线性近似贝叶斯估计相比其他估计的优越性,并考察了线性近似贝叶斯估计随着样本量增加的渐近性质.关键词:拉普拉斯分布;线性贝叶斯方法;Gibbs 采样;二次损失中图分类号:O212.8英文引用格式:YANG Y J,WANG L C.Approximate Bayesian estimation of the parameters ofLaplace distributionJ.Chinese J Appl Probab Statist,2024,40(1):1832.(in Chinese)1引言拉普拉斯分布通常又被称为拉普拉斯第一错误定律,由拉普拉斯于 1774 年提出,用于描述错误出现的频率和误差之间的函数关系.假设 X 服从参数为(,)的拉普拉斯分布,其密度函数为f(x;,)=12exp(|x|),其中 x R,R 为位置参数,0 为尺度参数.拉普拉斯分布的密度函数是不光滑的函数,其对数函数关于 求导,导数并不是处处存在的,所以拉普拉斯分布也不属于Cramer-Rao 正则族.这些数学性质都给拉普拉斯分布的参数估计带来了困难.一直以来,在统计文献中,拉普拉斯分布经常被用作其他(主要是正态)分布的反例,对该分布的研究也远远不及正态分布.然而,从图 1 中我们可以看到拉普拉斯分布具有尖峰厚尾的分布特点,所以当在数据中观察到比正态分布尾巴重的“东西”,它就会迅速成为首选分布.近年来,拉普拉斯分布被广泛地应用到许多领域,从图像、语音识别(输入分布)到海洋工程(导航误差的分布)和经济(商品对数收益的分布).Kotz 等1在The Laplace Distributionand Generalizations一书中收集整理了关于拉普拉斯分布理论以及应用研究的 400 多篇文献.国家自然科学基金项目(批准号:11371051)资助.通讯作者,E-mail:.本文 2021 年 1 月 4 日收到,2022 年 10 月 14 日收到修改稿.第 1 期杨彦娇,王立春:拉普拉斯分布参数的近似贝叶斯估计1910505100.00.10.20.30.40.5=0,=1=0,=2=0,=4=5,=4LaplaceNormal(0,1)图 1拉普拉斯分布和正态分布密度函数图参数估计方面,在大样本情况下,对拉普拉斯分布常用的有极大似然估计和矩估计.除此之外,Asrabadi2研究了位置参数已知时拉普拉斯分布函数的一致最小方差无偏估计,进而给出尺度参数的估计量;Rao 等3推导出了完备和截尾样本下拉普拉斯分布尺度参数的最佳线性无偏估计量.在小样本情况下通常使用贝叶斯方法,二次损失函数是贝叶斯方法下统计决策问题中用得最多的损失函数,记b 为参数 的估计,则其数学形式为L(b,)=(b )D(b ),其中 D 为正定矩阵.记 h(,|x)为(,)的联合后验密度,则在二次损失下,参数 和 的贝叶斯估计分别为b BE=h(,|x)dd,bBE=h(,|x)dd.因为 h(,|x)中存在绝对值的数学形式,这使得积分过于复杂,无法得到显示解.针对这种在贝叶斯估计中常出现的问题,一种处理方法是 Lindley4提出的近似计算两个积分比值的方法,它可以给出估计的近似数值解,但拉普拉斯分布并不满足其所需的正则条件;另一种常用方法是马尔可夫链蒙特卡罗方法(MCMC),但它往往需要迭代求解.因此,目前关于贝叶斯的相关研究主要局限在单参数或是没有显示解的估计上,如丁晓和韦来生5利用非参数方法构造了拉普拉斯分布位置参数的经验贝叶斯估计;徐美萍和段景辉6给出了单参数拉普拉斯分布的贝叶斯估计.和上述文献不同,本文利用线性贝叶斯方法构造了拉普拉斯分布的双参数的线性近似贝叶斯估计.该方法同时给出两个参数的估计,既20应用概率统计第 40 卷采纳了先验信息,但又不要求具体的先验分布,并且具有显示解,相对于迭代算法得到的贝叶斯估计具有稳健性的优势.本文安排如下:在第 2 节中,我们给出了拉普拉斯分布两参数的线性近似贝叶斯估计的表达式.第 3 节给出了线性近似贝叶斯估计相比其他估计的优越性的理论证明和数值模拟结果.第 4 节是结论部分.2线性近似贝叶斯估计线性贝叶斯方法最早由 Hartigan7提出,后由 Rao8从线性最优化角度进行了推广,随后线性贝叶斯方法又经过了很多贝叶斯统计学者的研究.如 Samaniego 和 Vestrup9使用线性经验贝叶斯方法来改进标准估计量;Wei 和 Zhang10将线性贝叶斯方法应用在线性回归模型中;周静雯和韦来生11将线性贝叶斯方法应用在生长曲线模型中;Zhang等12将对线性贝叶斯估计的研究扩展到多元线性模型;Wang 和 Singh13给出了型截尾下双参数指数族的线性贝叶斯估计.作为实际应用中贝叶斯估计的一种很好的近似,相比贝叶斯估计,线性近似贝叶斯估计相当于舍弃了全局最优,通过限定估计为线性结构,使得估计具有显示解.虽然线性贝叶斯方法是“局部最优”,但往往通过选择合适的统计量,可以得到相比 MCMC 方法更加稳定、便于使用的估计,因此统计量的选择至关重要.假设 X1,X2,Xn为来自拉普拉斯分布的简单随机样本.记 =(,2)是要估计的参数向量,T 为构造线性近似贝叶斯估计的统计量,且令T=(XS2/2),其中 X=(ni=1Xi)/n 和 S2=ni=1(XiX)2/(n1).易见有 E(X|)=,E(S2|)=22.定理 1T 的条件协方差阵为对角矩阵,即 Cov(X,S2|)=0.证明:Cov(X,S2|)=Covni=1Xin,ni=1(Xi X)2n 1?=1n(n 1)ni=1Cov(Xi,ni=1X2i?)1nCovXi,(ni=1Xi)2?.由于 E(Xi|)=,E(X2i|)=2+22和 E(X3i|)=62+3,因此Cov(Xi,ni=1X2i?)=Cov(Xi,X2i|)+nj=1,j=iCov(Xi,X2j|)第 1 期杨彦娇,王立春:拉普拉斯分布参数的近似贝叶斯估计21=E(X3i|)E(Xi|)E(X2i|)+(n 1)E(XiX2j|)E(Xi|)E(X2j|)=42.又有Cov(Xi,ni=1nj=i+1XiXj?)=(n 1)22,于是Cov(Xi,ni=1X2i?)1nCovXi,(ni=1Xi)2?=0,从而 X 和 S2不相关.?注记 2定理 1 说明了选择(X,S2)为T可以确保我们最终得到的估计不会过于复杂.而事实上,如果选择(X,S)为 T,则其条件协方差矩阵无法求得显示解.结合以上两点,我们设计估计参数为(,2)而不是(,).在实际应用中得到(,2)的估计后,再对第二项开根号就可以得到对 的估计.假设参数(,)的先验分布 G(,)满足 E(2+4)Cov(),从而得到MSEM(bM)MSEM(bLB)LWL+(L I)Cov()(L I)Cov()Cov(T)1Cov(T)Cov()=LCov(T)L 2LCov()+Cov()Cov(T)1Cov()=LCov(T)1/2 Cov()Cov(T)1/2LCov(T)1/2 Cov()Cov(T)1/2=UU 0.第 1 期杨彦娇,王立春:拉普拉斯分布参数的近似贝叶斯估计25定理得证.?因为均方误差矩阵准则下bLB优于bM是 MSE 准则下bLB优于bM的充分条件,所以我们可以得到结论:MSE 准则下bLB优于bM.数值模拟将在后面 3.3 节给出.3.2贝叶斯估计我们考虑以下情况的先验分布:面对只知道参数、的取值范围分别为 a,b,c,d 的情况时,我们可以取区间上的均匀分布作为其信息先验,这是由 Laplace14首先提出和使用的.则联合后验密度函数为H1(,|x)ne(ni=1|xi|)/,a b,c d.相应的贝叶斯估计为b BE=ne(ni=1|xi|)/ddne(ni=1|xi|)/dd,bBE=n+1e(ni=1|xi|)/ddne(ni=1|xi|)/dd.因为被积函数中有绝对值项,使得我们无法得到显示解.利用 Gibbs 取样方法,此时 和 的满条件分布为1(|,x)e(ni=1|xi|)/I(a b),2(|,x)ne(ni=1|xi|)/I(c d).通常情况下,计算没有显示解的贝叶斯估计常常先使用 Laplace 近似算法(参见文献15)获取初始值,再用该初始值进行 Gibbs 取样,这样可以减少 Gibbs 取样需要迭代的次数.因为对 f 取对数后 Hessian 矩阵不存在,导致这种情况下不能使用 Laplace 近似算法取初始值,所以我们使用无偏估计 T 作为 Gibbs 取样的初始值,进而得到贝叶斯估计.具体流程如下:1)给出两参数,的初始值,记为 0,0,并将第 j 步产生的 和 分别记为 j和 j;2)产生服从 1(|j,x)的 j+1;3)产生服从 2(|j+1,x)的 j+1;4)重复 2)和 3)N次;5)通过(N m0)1Nj=m0+1(j,j)计算(,)的贝叶斯估计值,其中 m0为无效样本的次数.26应用概率统计第 40 卷0.00.51.01.52.01.51.00.50.00.51.0log()图 2 和 ln()的联合后验密度的 0.1%、1%、10%的等高线及采样点图 2 显示 Gibbs 取样产生的参数样本点都落在了 10%的等高线内,即 MCMC 链是收敛的,因此 Gibbs 算法有效.由于贝叶斯估计没有显式解,我们通过数值实验来比较 MSE 准则下线性近似贝叶斯估计和贝叶斯估计的优劣,并考察线性近似贝叶斯估计随样本量增加的渐近性质.记eLB=(e LB,eLB),其中 e LB为bLB中的第一项,eLB为bLB第二项的算术平方根,bBE=(b BE,bBE)为(,)的贝叶斯估计.这里均匀先验下线性近似贝叶斯估计的表达式为bLB=n(b a)2X+4(a+b)(c2+cd+d2)8(c2+cd+d2)+n(b a)2n(n 1)(d c)2(4d2+4c2+7cd)S2+6(5n 3)(4i=0aib4i)(c2+cd+d2)2n(n 1)(d c)2(4d2+4c2+7cd)+18(5n 3)(4i=0aib4i).我们分别考察eLB和bBE的 MSE 随样本数和先验的分散程度变化而变化的情况.定义数值实验所用样本 MSE 为MSE=mi=1(b )(b )m,其中b 指的是估计值,为真实值,m 是估计的总次数.MSE 既考察了估计值和真实值的近似程度,又反映了估计本身的波动情况.实验一:考察随样本量 n 增加线性近似贝叶斯估计和贝叶斯估计的表现.取 a=1,b=1,c=0.5,d=1,则 和 分别来自 U(1,1)和 U(0.5,1),并作为 Laplace 分布的参数.样本量 n 的取值范围为 10 到 1000,间隔为 5.令 Gibbs 抽样迭代次数 N 为 3000,m0取 1000,m=50.第 1 期杨彦娇,王立春:拉普拉斯分布参数的近似贝叶斯估计27020040060080010000.000.050.100.15nMSEBayesLinear Bayesn(3/2)n(4/5)图 3实验一下eLB,bBE的 MSE 随样本量增加的变化从图 3 可以看出随着样本量的增加线性近似贝叶斯估计和贝叶斯估计的 MSE 十分接近,且均以介于 n4/5和 n3/2之间的速度逐渐趋于 0.表 1 展示了图 3 中 n=10,50,100,500,1000 的具体结果.表 1实验一下eLB、bBE的 MSEn10501005001000MSE(eLB)0.1164416870.0222325430.0147342310.0044918040.001634642MSE(bBE)0.1370493860.0223659310.0158156280.0025595080.0010747585实验二:考察当样本量 n=10,20,50,100 时,不同的先验分散程度下线性近似贝叶斯估计和贝叶斯估计的表现.考虑 来自 U(1,1),来自 U(c,d),取 c=0.5,d=0.6,0.65,9.05,10.其他参数设置同实验一.从图 4 可以看出随着先验的分散程度的增加,eLB和bBE的 MSE 都在增加.当 n=10,20 时,整体看来 MSE 准则下eLB的表现比较好;n=50,100 时,两者的表现好坏就无法从图中判别.鉴于图中所表现出来的这些特点,显然,只让 d 取固定几个值,来比较eLB和bBE的 MSE 从而得到结论是不恰当的.因此,我们对四种样本量情形下的估计的 MSE 分别取均值,记为 MMSE.表 2 详细展示了 n=10,20,50,100 时线性近似贝叶斯估计和贝叶斯估计的 MMSE,以及通过迭代算法得到bBE所花费的时间.结合实验一两种估计的表现,我们可以得出结论:在均匀先验下,当样本量较小时,线性近似贝叶斯估计的表现更好.相对地,样本量较大28应用概率统计第 40 卷0510152002468n=10MSEmethodBELB024602468n=20MSEmethodBELB0.00.51.01.52.002468n=50MSEmethodBELB0.000.250.500.751.0002468n=100MSEmethodBELB图 4实验二下eLB,bBE的 MSE 随先验的分散程度的变化表 2实验二下eLB,bBE的 MMSEnMMSE(eLB)MMSE(bBE)t(s)101.6962975.084915365.3201.1284921.748059396.7500.59037660.6310435391.71000.32138370.2642932404.0时,贝叶斯估计表现更好,但是值得注意的是,随着样本量增加,求解bBE所花费的时间也在增加.3.3数值实验最后,我们用数值实验展示在小样本、先验信息比较集中的情况下,线性近似贝叶斯估计eLB、贝叶斯估计bBE、极大似然估计bML和矩估计eM的估计效果.记(,)的矩估计、极大似然估计分别为eM=(X,ni=1(Xi X)22n),bML=(median(X),ni=1|(Xi median(X)|n).分别考虑以下两种情形:第一种仍为均匀先验,令 来自 U(1,2),来自 U(0.5,0.6);第 1 期杨彦娇,王立春:拉普拉斯分布参数的近似贝叶斯估计29第二种情形假设 的先验为正态分布,而 的先验为逆伽玛分布,即 N(0,0.1)、Inverse-Gamma(50,150).我们通过样本来估计(,),这个过程均重复 1000 次.02550750.00.20.4|LB|Frequency0501000.00.51.0|BE|Frequency02550750.00.20.40.60.8|ML|Frequency02550751000.000.250.500.751.00|M|Frequency图 5第一种先验下eLB、bBE、bML、eM与真实值的距离02550751001250.00.51.01.5|LB|Frequency02550751000.00.51.01.5|BE|Frequency0306090120024|ML|Frequency03060900246|M|Frequency图 6第二种先验下eLB、bBE、bML、eM与真实值的距离30应用概率统计第 40 卷图 5 和图 6 展示了四种估计与真实值距离的频数直方图.我们可以看到两种情形下线性近似贝叶斯估计离真实值的距离均明显小于极大似然估计和矩估计离真实值的距离.而与贝叶斯估计相比,线性近似贝叶斯估计与真实值的距离在 1000 次实验中有更多的频率落在 00.4 之间.因此,线性近似贝叶斯方法的估计效果最好.之前提到丁晓和韦来生5利用非参数方法构造了拉普拉斯分布位置参数的经验贝叶斯估计(BEB).进一步地,在图 7 中,我们模拟研究了针对位置参数 的线性近似贝叶斯估计 e LB和经验贝叶斯估计 b EBE的表现.0501000.00.51.01.5prior1:|LB|Frequency0501000.00.10.20.3prior2:|LB|Frequency01002003000.02.55.07.5prior1:|EBE|Frequency0501001502000246prior2:|EBE|Frequency图 7两种先验下 b EBE及 e LB与真实值的距离可以看到相比参数方法,非参数方法的经验贝叶斯估计效果并不是很好,这和非参数方法本身没有参数方法准确的情况是相符的.一般,非参数方法的经验贝叶斯估计在小样本的情况下较为不稳定,对于这类估计,我们通常研究其大样本性质.我们选取了 2021 年 2 月 8 日到 2021 年 7 月 9 日的上证指数收盘价,共 100 个数据,K-S 检验发现上证指数收益率服从拉普拉斯分布.因为收益率数据常常近似独立,我们每轮随机抽取其中 40 个做一次估计,每次估计中随机抽取 20 个用于估计参数的先验信息,剩下 20 个作为样本信息,一共进行 10 轮.在对参数的先验做估计时仍采用重复抽样的方法,从 20 个样本中随机抽样 500 次,每次抽取容量为 15 的一个样本,得到对 的估计.表 3 展示了在实际应用中线性近似贝叶斯估计和贝叶斯估计的表现情况.从均值来看,两种估计十分接近;从方差来看,线性近似贝叶斯估计要更稳定.这和图 4 中 n=20 时线第 1 期杨彦娇,王立春:拉普拉斯分布参数的近似贝叶斯估计31表 3eLB和bBE的表现e LBeLBb BEbBE均值0.0012600.0070860.000840.008546方差2.835751 1061.918982 1063.324537 1063.535534 106性近似贝叶斯估计的 MSE 小于贝叶斯估计的 MSE 的结论是一致的,也说明了有显示解的线性近似贝叶斯估计在小样本情况下的优越性.4结论本文针对传统贝叶斯方法下拉普拉斯分布两参数估计的困难提出了线性近似方法,给出其具有显式解的贝叶斯估计,并说明了选取 X、S2作为统计量的理由.文章给出了均方误差矩阵准则下,线性近似贝叶斯估计优于矩估计的理论推导,并用数值模拟的方法考察了均匀先验下线性近似贝叶斯估计和通常贝叶斯估计的 MSE 的表现:在小样本情况下线性近似贝叶斯估计优于贝叶斯估计,这意味着在应用中不需要迭代求解的线性近似贝叶斯更具有优势.在小样本且先验信息较为集中的情况下,线性近似贝叶斯估计与其他估计相比更接近真实值.此外,我们还发现随着样本量 n 逐步增加,线性近似贝叶斯估计的 MSE以不小于 n4/5的速度趋向 0.参考文献1 KOTZ S,KOZUBOWSKI T J,PODGORSKI K.The Laplace Distribution and Generalizations:ARevisit with Applications to Communications,Economics,Engineering,and FinanceM.Boston,MA:Birkh auser,2001.2 ASRABADI B R.The exact confidence interval for the scale parameter and the MVUE of the LaplacedistributionJ.Comm Statist ATheory Methods,1985,14(3):713733.3 RAO A V,RAO A V D,NARASIMHAM V L.Optimum linear unbiased estimation of the scaleparameter by absolute values of order statistics in the double exponential and the double WeibulldistributionsJ.Comm Statist Simulation Comput,1991,20(4):11391158.4 LINDLEY D V.Approximate Bayesian methodsJ.Trab Estad Investig Oper,1980,31(1):223245.5 丁晓,韦来生.双指数分布位置参数的经验 Bayes 估计问题J.数学杂志,2005,25(4):413420.6 徐美萍,段景辉.Laplace 分布参数估计的损失函数和风险函数的 Bayes 推断J.统计与决策,2010,(1):1314.7 HARTIGAN J A.Linear Bayesian methodsJ.J Roy Statist Soc Ser B,1969,31(3):446454.8 RAO C R.Linear Statistical Inference and Its ApplicationsM.2nd ed.New York:Wiley,1973.9 SAMANIEGO F J,VESTRUP E.On improving standard estimators via linear empirical BayesmethodsJ.Statist Probab Lett,1999,44(3):309318.32应用概率统计第 40 卷10 WEI L S,ZHANG W P.The superiorities of Bayes linear minimum risk estimation in linear modelJ.Comm Statist Theory Methods,2007,36(5):917926.11 周静雯,韦来生.生长曲线模型中参数的Bayes线性无偏估计J.应用概率统计,2008,24(6):639647.12 ZHANG W P,WEI L S,CHEN Y.The superiorities of Bayes linear unbiased estimator in multivariatelinear modelsJ.Acta Math Appl Sin Engl Ser,2012,28(2):383394.13 WANG L C,SINGH R S.Linear Bayes estimator for the two-parameter exponential family undertype II censoringJ.Comput Statist Data Anal,2014,71:633642.14 LAPLACE P S.Th eorie Analytique des Probabilit esM.Paris:Mme Ve Courcier,1812.15 ALBERT J.Bayesian Computation with RM.2nd ed.New York:Springer-Verlag,2009.Approximate Bayesian Estimation of the Parameters ofLaplace DistributionYANG YanjiaoWANG Lichun(School of Mathematics and Statistics,Beijing Jiaotong University,Beijing,100044,China)Abstract:The Laplacian distribution is one of the most important distributions used to characterizethe peak and thick-tailed data.This paper proposes a linear approximation Bayesian estimation withexplicit solutions for the two parameters of the Laplace distribution.The superiority of linear approximateBayesian estimation over other estimators is verified by theoretical derivation and numerical simulations,and the asymptotic behavior of the linear estimation with the increase of sample size is investigated.Keywords:Laplace distribution;linear Bayes method;Gibbs sampling;quadratic loss2020 Mathematics Subject Classification:62F15