维纳过程单变点模型的贝叶斯参数估计

维纳过程单变点模型的贝叶斯参数估计
何朝兵
【摘要】通过引入潜在变量，利用正态分布的重要性质得到了维纳过程单变点模
型比较简单的似然函数。

结合Metropolis-Hastings算法对参数进行Gibbs抽样，基于Gibbs样本对参数进行估计。

随机模拟的结果表明估计的精度较高。

%By introducing a latent variable, the simple likelihood function of Wiener process with a change-point is obtained according to the important property of the normal distribution. All the parameters are sampled by Gibbs sampler together with Metropolis-Hastings algorithm, and the parameters are estimated based on the Gibbs samples. Random simulation results show that the estimations are fairly accurate.
【期刊名称】《湖南师范大学自然科学学报》
【年(卷),期】2016(039)004
【总页数】5页(P84-88)
【关键词】潜在变量;可加性;满条件分布;Gibbs抽样;Metropolis-Hastings算法【作者】何朝兵
【作者单位】安阳师范学院数学与统计学院，中国安阳 455000
【正文语种】中文
【中图分类】O212.8;O212.4
变点问题成为近年来比较热的研究方向,它在经济、质量控制和医学等领域应用广泛[1-5].变点分析方法主要有非参数方法、最小二乘法和贝叶斯方法等.而随着统计计算技术的发展,贝叶斯变点分析方法越来越受到人们的欢迎,而复杂性的计算是贝叶斯方法的难点.贝叶斯计算方法中的Markov chain Monte Carlo (MCMC) 方法是最近发展起来的一种简单有效的计算方法.MCMC方法中的Gibbs抽样和Metropolis-Hastings算法使变点分析变得非常方便[6-9].Gibbs抽样可以简化变点问题,例如,未知参数的满条件分布可转化为无变点的后验分布,变点的满条件分布可转化为分布参数已知的后验分布.维纳过程是具有平稳独立增量的二阶矩过程,是一种特殊的扩散过程,它在纯数学、应用数学、经济学与物理学中都有重要应用.维纳过程不只是布朗运动的数学模型,在应用数学中,维纳过程可以描述高斯白噪声的积分形式;在电子工程中,维纳过程是建立噪音的数学模型的重要部分;控制论中,维纳过程可以用来表示不可知因素.对扩散过程变点模型的研究较多[10-13],虽然维纳过程是特殊的扩散过程,但对维纳过程变点模型的研究却较少[14-15],并且这些文献都是基于随机微分方程的求解来进行参数估计,计算比较繁琐,但基于似然函数并且利用MCMC方法研究此模型还不多见.
本文主要利用MCMC方法研究维纳过程单变点模型的参数估计问题.通过添加潜在变量得到了比较简单的似然函数,结合Metropolis-Hastings算法对参数进行Gibbs抽样,基于Gibbs样本对参数进行估计.随机模拟的结果表明估计的精度较高. 定义1 随机过程W(t)如果满足:
1)W(0)=0,具有独立增量;
2)对任意s,t>0,W(s+t)-W(s)服从正态分布N(0,σ2t),σ>0,则称W(t)为以σ2为参数的维纳过程.
当维纳过程的参数σ2在某个时刻改变时,有如下定义.
定义2 随机过程W(t)如果满足:
1)W(0)=0,具有独立增量,
2)对任意
,则称此模型为维纳过程单变点模型,τ称为变点位置参数.
在n个时刻t1<t2<…<tn观察维纳过程单变点模型,得到观察数据
W(ti)-W(ti-1)=zi,t0=0,i=1,2,…,n.
假设已知在观察时间区域(0,tn]内有一个变点,即0<τ<tn,则必存在m使
τ∈(tm,tm+1],0≤m≤n-1,实际上m是τ的函数.
由式(1)知
).
由正态分布的可加性得
).
W(tm+1)-W(tm)的观察值为zm+1.
关于参数似然函数比较复杂,为了得到较简单的似然函数,不妨引入潜在变量
X=W(τ)-W(tm),下面求在W(tm+1)-W(tm)=zm+1的条件下,X的条件密度函数. 下面介绍概率论中一个很重要的结果,即下面的引理1.
引理1 设独立,令Z=X+Y,则在Z=z条件下,
,其中.
由引理1知,,
其中.
令Z,T分别表示由zi,ti组成的向量.记{m+2,…,n}.
当m=0时,D1=∅, 当m=n-1时,D2=∅.记x为X的取值,添加潜在变量后的似然函数为
,
其中∅时,sj=0,j=1,2.
下面给出参数的先验分布.
1) 取τ的先验分布为均匀分布,即π(τ)∝1,0<τ<tn.
2) 取的先验分布为倒伽玛分布IGa(αi,βi),即
.
假设,τ相互独立,则
(τ).
而x的满条件分布则由正态分布抽样,即).
简记的满条件分布为|·).,τ的满条件分布如下:
,
,
g(τ),0<τ<tn.
x的满条件分布为正态分布,的满条件分布为倒伽玛分布,所以都可以利用统计软件直接进行Gibbs抽样.但τ的满条件分布比较复杂,可以利用 MCMC方法中的Metropolis-Hastings算法对其抽样,此时选择均匀分布作为建议分布.
下面介绍MCMC方法的具体步骤.
设参数的初始值),第t-1次迭代结束时的估计值为θ(t-1),则第t次迭代分为如下几步:
1)令θ=θ(t-1),根据τ确定m,由随机产生x;
2)由抽取,用更新
3)由抽取,用更新
4)首先根据步骤1)中的m,x计算g(τ(t-1)),然后由均匀分布U(0,tn)抽取变点τ′,根据τ=τ′确定新的m,由新的随机产生x,根据新的m,x计算g(τ′).令,若U(0,1)随机数u≤α(τ(t-1),τ′),则τ(t)=τ′,否则τ(t)=τ(t-1),用τ(t)更新τ.
)称为θ的一个Gibbs样本.假设进行M次Gibbs抽样,第B次迭代以后抽样收敛,
根据后M-B个迭代值对参数进行估计.
下面进行随机模拟试验.
取(3,0.5,1.5).为了方便计算,等间隔进行观察W(t),即Δti=tn/n=Δ,此时
ti=iΔ,i=1,2,…,n.经计算得m=37,由分布产生m个随机数z1,…,zm,由分布产生1个随机数zm+1,由分布产生n-(m+1)个随机数zm+2,…,zn,则z1,z2,…,zn即是模拟的n个观察数据.根据这n个数据估计.
取的先验分别为IGa(5.5,2.5),IGa(2.3,3.5),取B=10 000,M=20 000.参数估计结果见表1.我们重点估计分析参数τ.τ的Gibbs抽样迭代过程见图1.Gibbs抽样收敛性诊断最常用的方法是产生多条马氏链,如果这几条链最后趋于重合,则抽样收敛.我们产生2条马氏链,τ的2条迭代链见图2.
由表1可以看出,把Gibbs样本的均值作为参数的估计时,τ的相对误差小于的相对误差不超过6%,估计精度较高,MC误差都很小. Gibbs样本的中位数与均值差别很小,所以把样本中位数作为参数的贝叶斯估计的精度也较高.各参数可信水平为0.95的可信区间可近似取为[2.5%分位数, 97.5%分位数],可以看出近似可信区间的长度非常短,所以区间估计的效果也较好;由图1可以看出迭代的波动很小.由图2可以看出τ的两条迭代链都趋于重合,说明抽样收敛.综上分析,随机模拟试验的效果较好.
【相关文献】
[1] PAGE E S. Continuous inspection schemes[J]. Biometrika, 1954,41(1):100-115.
[2] CHERNOFF H, ZACKS S. Estimating the current mean of a normal distribution which is subjected to changes in time[J]. Ann Math Stat, 1964,35(3):999-1018.
[3] CSÖRGÖ M, HORVTH L. Limit theorems in change-point analysis[M]. New York: Wiley, 1997.
[4] PERREAULT L, BERNIER J, BOBÉE B, et al. Bayesian change-point analysis in hydrometeorological time series. Part 1. The normal model revisited[J]. J Hydrol,
2000,235(3):221-241.
[5] FEARNHEAD P. Exact and efficient Bayesian inference for multiple changepoint problems[J]. Stat Comput, 2006,16(2):203-213.
[6] LIANG F, WONG W H. Real-parameter evolutionary Monte Carlo with applications to Bayesian mixture models[J]. J Am Stat Assoc, 2001,96(454):653-666.
[7] LAVIELLE M, LEBARBIER E. An application of MCMC methods for the multiple change-points problem[J]. Sig Pro, 2001,81(1):39-53.
[8] KIM J, CHEON S. Bayesian multiple change-point estimation with annealing stochastic approximation Monte Carlo[J]. Comput Stat, 2010,25(2):215-239.
[9] YUAN T, KUO Y. Bayesian analysis of hazard rate, change point, and cost-optimal burn-in time for electronic devices[J]. IEEE Trans Rel, 2010,59(1):132-138.
[10] ABBAS-TURKI L A, KARATZAS I, LI Q. Impulse control of a diffusion with a change point[J]. Stoch Int J Probab Stoch Process, 2015,87(3):382-408.
[11] MISHRA M N, PRAKASA RAO B L S. Estimation of change point for switching fractional diffusion processes[J]. Stoch Int J Probab Stoch Process, 2014,86(3):429-449. [12] GAPEEV P V, SHIRYAEV A N. Bayesian quickest detection problems for some diffusion processes[J]. Adv Appl Probab, 2013,45(1):164-185.
[13] NEGRI I, NISHIYAMA Y. Asymptotically distribution free test for parameter change in
a diffusion process model[J]. Ann Inst Stat Math, 2012,64(5):911-918.
[14] VOSTRIKOVA, L YU. Detection of a “disorder” in a Wiener pro cess[J]. Theor Probab Appl, 1981,26(2):356-362.
[15] HADJILIADIS O, MOUSTAKIDES V. Optimal and asymptotically optimal CUSUM rules for change point detection in the Brownian motion model with multiple alternatives[J]. Theor Probab Appl, 2006,50(1):75-85.。