MCMC

2024-05-30

MCMC（通用8篇）

MCMC 篇1

0 引言

通信信号调制方式的识别是软件无线电一个重要的方面, 已经广泛应用于电子对抗、情报截获等军事通信以及通信监视、频谱管理等民用通信方面。调制识别的基本方法, 文献将其分为两类:统计模式识别方法和似然比检验方法。统计模式识别方法是通过提取能反映信号调制类型的统计特征来进行模式识别, 常用的分类特征包括:信号的瞬时幅度、频率及相位信息, 高阶矩及高阶累积量, 功率谱及高阶谱特征等。在实际工程中, 要选择一个分类意义上的好的特征很困难。在大多情况下, 特征的选择依赖于系统设计者的直觉。同时, 要严格证明这些算法在统计意义上的最优也很困难。

最大似然理论是解决调制识别问题的另一个经典方法, 它在加性高斯白噪声环境下, 将调制识别问题转化为复合假设检验问题。由检测理论可知, 给定观测数据和备择的调制类型集合, 在贝叶斯最小误判概率准则下, 最优的调制识别分类器是通过最大化后验概率来实现的。当备择的各个调制类型先验等概率时, 最优分类器简化为最大似然分类器, 通过对信号的似然函数进行处理, 从备择假设中选择具有最大似然概率者为真。在电子战等非协作通信环境中, 由于数字通信的信息内容未知, 以及存在信号参数的估计误差, 构造的信号似然函数中一般含有未知参数, 此时, 最大似然统计假设检验通常采用平均似然比检测 (ALRT) 方法, 但是在计算过程中无法避免的多重积分问题使得该类算法的应用受限。

为保持ALRT算法的优越性, 同时解决上述多重积分计算复杂甚至大多情况下不存在解析解的问题, 文献采用Metropolis-Hastings (MH) 算法对ALRT中的多重积分进行数值估计。MH算法是一种基于数值计算的马尔可夫链蒙特卡罗 (MCMC) 方法, 在统计估计计算中越来越受到关注, 但是该方法中的建议分布函数的选取具有任意性, 如何选择最优的建议分布函数至今仍是难题。本文引入一种自适应的MCMC方法——Adaptive Metropolis (AM) 方法, 其建议分布函数在迭代过程中自适应的更新, 由此避免了建议分布函数的选取问题, 并且可以得到比MH算法更优的估计效果。

1 似然比方法

假设接收信号的复基带形式如下

其中备择调制类型有K种, {an (k) }为选自第k个调制星座A (k) ={a (k, 1) , (42) , }, k (28) 1, (42) , K的码元序列, Mk为该星座中的星座点个数, 并且不是一般性设。假设未知的频偏和相偏参数在码元观测时间内保持不变, 分别记为和f0。N为观测码元序列长度, {Vn}是方差为σ2的高斯白噪声序列。

基于似然函数的调制识别问题旨在通过对长度为N的被噪声污染的接收码元序列进行观测, 从K种可能的调制样式中选择出正确的调制类型, 它可被视为一个多元假设检验问题:

其中r (28) [r1, (42) , rN]T, 在噪声建模为高斯分布的假设下, 第k种假设的条件似然函数可表示为:

其中为由未知频偏和相偏组成的未知参数向量。平均似然比检测 (ALRT) 将未知参数看成概率密度函数 (PDF) 已知或者可以假设的随机变量, 似然函数是在对未知参数取平均意义上的结果, 可表示为:

其中为Hk假设下未知参数向量u的先验概率密度函数。当所有调制样式先验等概时, ALRT理论可通过选取具有最大似然概率的备择假设提供最佳识别效果。然而上式中的多重积分往往不存在闭式解, 另外, 多重积分的计算量随着向量u维数的增多而呈指数增长, 这些因素使得ALRT算法的应用受到极大限制。式 (3) 的求解问题是似然分类的关键问题, 本文采用基于数值计算的马尔科夫链蒙特卡罗 (MCMC) 方法使之得以有效实现。

2 基于AM算法的调制识别原理

近年来, 起源于统计物理学的MCMC方法引起了信号处理领域的极大关注, 其核心思想是产生一个以目标概率密度函数为稳定分布的各态历经的马尔科夫链, 并以该链达到稳态时的样本作为目标分布的样本。

根据贝叶斯理论, 很容易得到:

由重要采样理论, 式 (3) 的估计可以由下式得到:

其中ui是对目标分布函数 (uHk, r) 取样得到, Ns为用于式 (3) 的近似估计的取样点的个数。由于未知参数向量u的先验概率不依赖于假设Hk, 因此概率函数可简化为: (uiHk) (28)  (u i) 。取样点的遍历性保证了在取样个数Ns的情况下, 式 (5) 的数值无限接近于p (rHk) 。因此产生服从目标分布 (uHk, r) 的取样点是MCMC分类器的关键。

Metropolis-Hastings (M-H) 算法可产生服从目标分布的马尔科夫序列u0, u1, u2, (42) 。假设当前状态ui (28) u, 下一状态通过从建议分布函数q (|u) 中抽取待选采样点w产生。然而建议分布函数的选取以及调整方向很难确定, 为避免这一难题, 本文采用自适应Metropolis (AM) 算法。

AM算法是Haario在2001年提出的一种改进MCMC抽样算法。与传统M-H算法相比, AM算法不需要预先确定参数的建议分布函数, 而是根据目前已得到的所有取样点信息对建议分布函数进行自适应更新。

假设已抽取采样点u0, u 1, (42) ui, 下一个待选采样点w由建议分布函数qi (|u0, u 1, (42) ui) 中抽取产生, 并且w的接受概率为:

也就是说, 下一个取样点以概率a (u i, w) 被设置为ui (10) 1 (28) w, 以概率1-a (u i, w) 保持前一取样点的值不变, 即ui (10) 1 (28) ui。AM算法中的建议分布函数就是均值为为前i个取样点的均值向量) , 协方差矩阵为Ci (28) Ci (u 0, u 1, (42) ui) 的高斯分布函数。根据大数定义理论, 多个概率密度分布未知的参数的联合概率密度定义为高斯分布是较为合理的。选择合适的过渡迭代次数N0 (29) 0, 则协方差矩阵定义为:

其中参数sd仅取决于未知参数向量ui的维数d,  (29) 0是一个非常小的数, 它的存在是为了保证iC在迭代过程中的非奇异性。Id为d维单位阵。为减少计算协方差矩阵过程中的计算开销, 可以方便的计算iC的值, 有下面的递归表达式:

其中。文献证明了经过很短的过渡迭代后, 由AM算法得到的马尔科夫序列能很好的收敛于目标分布, 并且AM算法中的自适应动态更新的高斯建议分布函数比传统M-H算法中预先指定的固定不变的建议分布函数更能得到满意的效果。

AM算法的计算流程如下:

Step 0. (初始化) 任意选取初始状态u0及初始协方差矩阵C0, 设置i=0。

Step 1.采用下列步骤生成ui+1。

(1) 根据式 (8) 计算Ci。

(2) 产生服从正态分布的待选取样点

(3) 产生服从均匀分布的随机数, 根据式 (6) 计算待选取样点w的接受概率。若v≤ (u i, w) , 则接受ui (10) 1=w, 否则ui+1=ui。

Step 2.重复Step 1直到产生足够多的取样值可以满足计算式 (5) 的要求。

Step 3.对于每一种可能的调制样式假设Hk, k=1, (42) , K重复step 0~2。比较每一种假设下的似然函数值, 判定使似然函数达到最大者为真。

3 仿真结果

为验证AM分类器的性能, 这里做了如下仿真试验。仿真中假设备择调制样式包括8PSK, 16QAM, 4PAM和BPSK, 且各种调制样式先验等概。信号经过匹配滤波、载波恢复及码元同步等处理后, 其基带复包络采样序列形式如式 (1) 所示, 并且信号功率经过归一化处理, 即Es (28) 1。假设频偏和相偏统计独立且各自服从定义域内的均匀分布

于是式 (5) 中的先验概率分布函数p (u|Hk) 可表示为:

AM采样器中参数设置如下:

(1) 参数

(2) 过渡迭代次数:Nb=400;

(3) 稳态后用于计算式 (5) 中遍历均值的取样长度:

(4) 马尔可夫链的初始状态:

3.1 仿真1算法的收敛性仿真

取样点收敛于平稳分布的速度以及收敛程度是MCMC算法的重要考虑因素。该仿真对AM算法的收敛性进行了验证, 假设接收序列为BPSK信号, 信噪比为5d B。码元观测长度为N=100, 频偏和相偏参数分别设置为f0=0.2 N, θ0=p12。利用AM算法迭代产生的θ0和f0的马尔可夫采样序列如图1所示。由于初始状态的任意性, 马尔可夫链需要经过一定次数的迭代后才能达到稳定, 因此初始阶段取样值与参数真实值相差较大。随着迭代过程的演进, 采样点逐渐逼近参数真实值, 大约400次迭代后, 采样值几乎与真实值相同。相比于文献中的M-H算法, 本文算法收敛速度稍慢, 这是由于文献算法将未知参数进行分块处理, 每次迭代过程分别针对θ0和f0进行更新, 而本文算法将未知参数作为整体向量更新, 在迭代过程中θ0和f0的收敛相互受到制约, 因此影响了该算法的收敛速度, 然而相比于M-H算法, 该算法对参数估计的准确程度较高。

(θ0和f0的真实值分别为θ0 (28) 0.2617, f0=2×10-3)

3.2 仿真2算法分类性能比较

这里将本文算法与另外两种算法进行比较:文献中的M-H算法和最常见的基于4阶累积量的HOS算法。

图2给出了在SNR=5d B, N=100条件下, 载波频偏分别对三种算法分类性能的影响, 其中f0的变化范围为0~0.5 N。可以看出存在载波偏移的情况下, 本文所提出的AM分类器性能明显优于M-H分类器和HOS分类器。这是由于在经过短暂的过渡迭代后, 频偏和相偏都非常接近于其真实值, 如图1所示, 因此算法对频偏和相偏具有鲁棒性。HOS分类器的正确识别概率随载波频偏的增大急剧下降, 其性能受频偏的影响最大, 而AM算法对频偏的鲁棒性最强。

由于HOS算法在提出时并未考虑频偏和相偏问题, 它仅仅适用于零频偏和零相偏的较理想情况, 于是这里假设在进行HOS计算之前的预处理阶段已经通过某些算法实现了频率和相位的精确估计。为验证本文所提出的AM算法的优越性, 这里将频偏和相偏分别为θ0=p12, f0=0.2 N的AM算法和M-H算法与零频偏、零相偏的HOS算法进行比较。图3为各算法分别对8PSK/16QAM/4PAM/BPSK四种信号进行1000次调制识别试验的平均正确识别概率 (PCC) 随信噪比 (SNR) 变化曲线。试验中信噪比的变化范围是0d B到10d B, 变化步长为1d B, 码元个数N=100。仿真结果表明即使在不存在频偏与相偏的情况下, HOS分类器的识别效果在三种分类器中表现最差, 相比M-H算法, AM算法的正确识别概率有所提高, 在SNR为5d B时, 正确识别概率大于95%, 验证了算法的有效性。

3.3 仿真3运算量比较

本文提出的AM算法与传统M-H算法的运算量分析如下。由于两种算法在每一次迭代中均需生成正态分布及均匀分布的随机数, 此过程可由多种实现方法, 例如可以采用逼近法、反函数法、舍选法等生成正态分布随机数。不同随机数生成算法将导致不同的运算量, 因此这里只能对AM与M-H算法的运算量做定性分析。由于M-H算法将未知参数分块化处理, 在本例中0和f0两个参数在每次迭代中单独更新, 即每一次迭代需分别完成两次提取服从正态分布的取样点与计算接收概率的过程, 而AM算法将未知参数进行整体更新, 在每一次迭代中仅需完成一次该运算, 同时, 虽然建议分布函数在每一次迭代中需实时更新, 由于式 (8) 递归形式的引入, 使这一运算量相对很小。综上所述, AM算法的总体运算量小于M-H算法。

4 结论

本文提出了一种基于似然函数的调制识别算法, 利用自适应Metropolis (AM) 算法, 产生满足目标分布的未知参数和发送符号的各态历经样本从而实现似然函数的近似计算。相比于传统Metropolis-Hastings算法, 该算法避免了建议分布函数选取这一难题。仿真结果表明本文算法识别效果优于Metropolis-Hastings算法以及高阶累量算法, 证明了该算法的有效性。

参考文献

[1]O.A.Dobre, A.Abdi, Y.Bar-Ness and W.Su.Survey of automatic modulation classification techniques:classical approaches and new trends[J].IET Commun.2007.

[2]Anna, K and David, K.Algorithms of digital modulation classification and their verification[J].WSEAS Transactions on Communications.2010.

[3]A.Swami and B.Sadler.Hierarchical digital modulation classification using cumulants[J].IEEE Trans.Comm.48 (3) .2000.

[4]Zhang L.P and Wang J.X.Blind modulation recognition algorithm for MQAM signals.Journal of Electronics&Information Technology.2011.

[5]D.C.Popescu, O.A.Dobre and F.Hameed.On the likelihood-based approach to modulation classification[J].IEEE Trans.Wireless Communications.2009.

[6]S.Lesage, J.Y.Tourneret and P.M.Djuric.Classification of digital modulations by MCMC sampling[J].IEEE ICASSP.2001.

[7]柳征, 王明阳, 姜文利等.一种新的贝叶斯调制分类算法[J].电子与信息学报.2006.

[8]S.Chib and E.Greengerg.Understanding the Metropolis-Hastings algorithm.American Statistician.1995.

[9]H.Haario, E.Saksman and J.Tamminem.An adaptive Metropolis algorithm.Bernoulli.2001.

MCMC 篇2

关键词马尔科夫链蒙特卡洛模拟；贝叶斯估计；改进广义帕累托分布；地质灾害

中图分类号 O213.2 文献标识码 A

The MGPD Model Based on MCMC Simulation

and Its Application in Geological Disaster Risk Measure

OUYANG Difei1，YANG Yang1， GAN Liu2， LI Yingqiu1

（1.School of Mathematics and computing Science， Changsha University of Science

and Technology， Changsha，Hunan 410004，China；

2. School of treasury and finance， Hunan University of Commerce， Changsha， Hunan 410205，China）

Abstract We used Bayesian estimation based on Markov Chain simulation to optimize the meliorated Generalized Pareto Distribution model （MGPD）， and obtained the estimation of the Value at Risk（VaR） and Conditionl Value at Risk（CVaR）. The empirical study and adaptability test of the model were based on geological disasters loss data of Loudi City in Hunan Province. The conclusion shows the optimized model has not only excellent ability in describing the data， but also extensive applicability.

Key words Markov Chain Monte Carlo simulation； Bayesian estimation； meliorated generalized Pareto distribution model； geological disaster

1 引言

地质灾害是指在地质作用下，地质自然环境恶化，造成人类生命财产损毁或人类赖以生存与发展的资源和环境发生严重破坏的过程或现象.据国土资源部统计，2013年，全国共发生各类地质灾害15403起，造成481人死亡、188人失踪、264人受伤，造成直接经济损失101.5亿元.死亡人数与上年相比，同比增加7.5%.地质灾害风险评估作为一项极具现实意义的重要研究课题和减轻灾害损失的非工程性重要措施，其研究成果已经引起了社会的广泛关注.这其中涉及一系列与统计理论相关的方法，通过对地质灾害风险进行评估及管理来刻画地质灾害风险，对政府及保险机构防范风险、稳定经营、降低破产概率就显得至关重要.

地质灾害风险导致索赔的统计数据数量不多、质量不高，因此在进行风险研究时采用传统的精算方法很难准确预测未来损失和管理风险.极值理论常用于研究随机变量，或一个随机过程的随机性质，最常见的是指在特殊情况下发生极端事件的概率.因此，在分析解决地质灾害风险等随机问题时，极值理论大有用武之地.极值统计中主要有两类模型，一类是区块极值模型（BlockMmaximum Model，简称为BMM），这种模型主要对组最大值建模.另一类是基于广义帕累托（Generalized Pareto）分布的模型（简称GPD模型），它是对观察值中所有超过某一阈值的数据建模.由于GPD模型充分利用了数据中的极值信息，因此针对地质灾害风险导致统计数据数量不多、质量不高的情况，采用GPD模型将更有用[1].

经济数学第 32卷第2期

欧阳迪飞等：基于MCMC模拟的MGPD模型及其在地质灾害风险度量中的应用

以湖南省娄底市地质灾害为例，利用MGPD实证得到了当地的地质灾害损失的在险风险值和条件在险风险值.首先，根据QQ图和经验平均超出函数图对原始数据进行诊断并选取阈值，结果显示样本数据具有厚尾的特征；其次，选择扩展Burr XII分布对这类数据进行描述，并基于MCMC贝叶斯估计确定分布函数中的参数估计值；最后，利用所得到的分布函数检验了模型的适应性以及测算了不同置信水平下的在险风险价值损失、条件在险风险价值损失，并据此说明了研究结论和实际意义.

2 MGPD模型

2.1 扩展Burr XII分布

MGPD模型将所有超出给定充分大阀值的观测值作为观测样本，进而研究观测值的渐进分布.MGPD模型基于扩展Burr XII分布.

4 结论

首先，MGPD模型能很好地描述地质灾害损失数据，对尾部极值数据的捕获能力较高.其次，基于MCMC模拟的贝叶斯估计对模型的适应性具有较好地改善，能够保证一定程度下的数据波动不会对模型造成明显的干扰.其次，实证得出在99%置信水平下认为未来某一次的地质灾害在险风险损失为197.103万元，显然，为规避巨灾损失，我国势必需要加强对地质灾害的防治与预警力度.

nlc202309031229

地质灾害在中国乃至全世界都是一个无法回避的问题，如何有效地规避地质灾害风险，尽可能的降低地质灾害给国家、政府以及人民带来的影响是一个在很长时间内都需要面对的问题.一方面应该做好地质灾害的预防工作，避免因为相关设施的落后、反应机制的不健全，造成不必要的人员财产损失.另一方面，即在灾害发生之后，如何有效、高效地减轻地质灾害带来的负面影响，于国于民无疑是有重要意义的[10].显然，只有综合考虑这两个方面，才能防患于未然，我国有必要加快地质灾害防治体系的完善，尽可能降低地质灾害给国民经济和人民生活带来的不利影响.

参考文献

[1] 欧阳资生. 地质灾害损失分布拟合与风险度量[J]. 统计研究， 2012， 28（11）： 78-83.

[2] Shao Q， Wong H， Xia J， et al. Models for extremes using the extended threeparameter Burr XII system with application to flood frequency analysis[J]. Hydrological Sciences Journal， 2004， 49（4）： 685-702.

[3] 欧阳资生. 极值估计在金融保险中的应用[M]. 北京：中国经济出版社， 2006： 127-129.

[4] I USTA. Different estimation methods for the parameters of the extended Burr XII distribution[J]. Journal of Applied Statistics， 2013， 40（2）： 397-414.

[5] G MATTHYS， J BEIRLANT. Estimating the extreme value index and high quantiles with exponential regression models[J]. Statistica Sinica， 2003， 13（3）： 853-880.

[6] 刘睿，詹原瑞，刘家鹏. 基于贝叶斯MCMC的POT模型——低频高损的操作风险度量[J]. 管理科学， 2007， 20（3）： 76-83.

[7] 李应求，田琴，戴志锋. 基于鲁棒均值下半偏差模型的供电公司购电组合策略[J]. 经济数学， 2014， 31（4）： 14-19.

[8] 李应求，甘柳，魏民. 一类多险种复合PoissonGeometric过程风险模型研究[J]. 统计与决策， 2010 （7）： 53-55.

[9] 李应求，刘薇，陈文锋. 聚类分析视角下地区保险业发展差异研究—基于湖南省各地市的截面数据分析[J]. 时代金融， 2009（1）：117-119.

[10]李应求，刘朝才，彭朝晖. 不确定条件下企业的投资规模决策[J]. 运筹学学报， 2008， 12（2）： 121-128.

MCMC 篇3

关键词：大规模MIMO,MCMC检测,超松弛迭代,Gibbs采样

0 引言

大规模MIMO(Large-scale MIMO) 又称Massive MIMO,最早由美国贝尔实验室的MARZETTA T L提出[1], 该技术在基站配置数十根甚至上百根天线,以获得更大的空间自由度。文献[1-4]研究表明,当天线数目趋于无穷大时,瑞利衰落和加性高斯白噪声等负面影响可以忽略不计,数据传输速率能得到极大提高。大规模天线阵列既带来了性能增益, 也带来了前所未有的挑战, 如传输方案设计、迅速增加的硬件复杂度和计算量等。大规模MIMO系统中低复杂度、有效的接收算法是该技术从理论到实现的关键。

近年来, 统计学中的一些方法, 如基于蒙特卡罗仿真的MCMC检测算法[5,6,7,8]被引入到MIMO系统中, 它利用已知符号的后验条件概率来迭代地求出下一个符号的后验条件概率, 最后得到发送矢量的后验概率。MCMC方法的性能主要取决于采样点数量和迭代次数,与待估计变量维数无关,可以避免算法复杂度随天线数和调制阶数呈指数级增长的问题。传统的MCMC算法需要经过足够次数的采样后,马尔可夫链才能趋于平衡分布。另外,在高信噪比条件下容易出现“陷入”问题(stalling problem )[7], 即采样 “ 陷入” 某一固定状态, 采样状态减少,导致后验概率估计误差。针对上述问题,本文采用超松弛迭代的方法构造马尔可夫链,选择合适的松弛因子加快马氏链的收敛速度。仿真结果表明,该算法提高了系统检测性能,降低了运算复杂度。

1 系统模型

本文研究的大规模MIMO系统由一个天线数为N的基站和K个单天线用户构成,K≤N,考虑该系统的上行链路,基站端的接收向量表示为:

其中x为K个用户的信号发送向量,x=[x1, x2, … , xK]T;H为N × K维信道矩阵; n是均值为零、协方差N0IN的加性高斯白噪声,n=[n1, n2, … , nN]T。

最大似然(Maximum Likelihood ,ML) 检测是MIMO系统的最优检测算法, 通过遍历所有可能的发送信号组合,寻求最优的检测值,即:

式中(χ)K表示发送向量x的全部可能取值。随着收发天线数及调制阶数的增加,ML算法搜索空间呈指数级增长,实际难以实现。

2 大规模MIMO信号检测算法

大规模MIMO基站天线数量可能达到几百根以上,传统MIMO的一些准最大似然方法( 如基于QR分解的QRM - MLD算法和球形译码( SD ) 算法) 在这里难以采用。探索大规模MIMO系统中高性能、低复杂度的信号检测算法是要解决的主要问题。文献[9] 将MCMC方法引入大规模MIMO,并获得了较好的性能。

2 . 1 MCMC算法

MCMC检测通过统计抽样获得发送符号矢量, 用统计方法估计各符号的后验概率,其性能和运算量与发送信号的维数无关,只取决于采样点数量和迭代次数。

MIMO系统中, 多维信号的联合概率分布如下:

MCMC算法从条件分布p ( x | y , H ) 中抽取样值x , 形成马尔可夫链。 MIMO系统中马尔可夫链状态数随x的维数呈指数增加, 为了降低采样复杂度, 采用Gibbs采样构造马尔可夫链。 Gibbs采样[10]是一种基于条件分布的迭代采样方法,它利用已知符号的后验条件概率迭代地求出下一符号的后验条件概率,最后得到发送信号矢量的后验概率。第t+1 次迭代中第k个符号的Gibbs采样过程如下:

初始采样序列x(0)随机产生。由于:

令,Gibbs采样可以写为:

其中Ck是归一化参数,χ 是全部发送信号集合。

2 . 2 MCMC改进算法

传统MCMC算法需要经过足够多次迭代才能收敛至平衡分布,提高马尔可夫链收敛速度能够降低检测算法的计算复杂度[11], 这对大规模MIMO系统尤为重要。因此考虑用超松弛迭代方法(Successive Over-Relaxation ,SOR )[12]来提高马尔可夫链收敛速度。

将MIMO迭代检测器看作一个随机线性系统, 考虑基本线性系统模型如下式:

其中,M=HHH,b=HHy。当K取值很大时,M矩阵求逆运算复杂度很高,考虑采用迭代方法来求解[13], 由于M是对称矩阵,有:

其中,D是对角矩阵,L是严格下三角矩阵。由Jacobi迭代得到第t+1 次迭代的解x(t+1):

超松弛迭代法(Successive Over-Relaxation ,SOR) 引入了松弛因子 ω,以加快迭代矩阵的收敛速度,令:

当 ω∈(1,2)时为超松弛。 SOR迭代可以写作:

其中,。天线k发送信号更新为:

其中,

由系统模型(1)得:

将上式与式(11)比较,有:

其中,w=-HHn,是均值为零、协方差N0HHH的加性高斯白噪声。因此可以由以下分布得到采样值:

其中,Ck是归一化参数,σk2= N0| | hk| |2, 且:

SOR - MCMC检测算法实现过程如下:

(1)输入接收向量y、信道H、迭代次数T;

(2)随机产生初始向量x(0);

(3)while t<T do

(4)for k=1:K

( 5 ) { 由式( 16 ) 、 ( 15 ) 计算由SOR迭代获得的采样点x(t+1)k}

( 6 ) 由获得的采样值, 进行对数似然比( LLR ) 计算,并进行软判决。

2 . 3 算法复杂度分析

由上面的分析可知,MCMC算法的复杂度主要集中在由迭代采样构造马氏链,可以先不考虑预条件矩阵处理,即计算和b的运算量。由式(6)可知,传统MCMC算法获得采样值的复杂度为, 完成一次迭代的复杂度为; 而由式(15) ,SOR -MCMC中获得采样值的复杂度仅为, 完成一次迭代复杂度为。当迭代次数增加时, SOR - MCMC复杂度比MCMC算法低更多, 特别适用于大规模MIMO系统。

3 仿真结果及分析

接下来分析上述算法在不同天线方案、不同调制方式下的检测性能。大规模MIMO系统上行链路收发天线数分别为N和K ,记为K×N ,令收发天线比, 要求K ≤ N 。为了便于比较算法性能, 采用单用户匹配滤波器检测(MFD) 算法近似最佳的MLD算法, 即只考虑单个用户发送数据,利用MF方法恢复发送信号。

首先考虑SOR-MCMC算法中松弛因子 ω 对误码率(BER) 性能的影响。采用4QAM调制,K =16 , 信噪比SNR=10 d B时, 迭代次数T取20 , 基站天线数分别为8 / 16 / 32的仿真结果如图1。可以看到,β 较小时,ω 对误码率的影响更明显。 β 较大时,由于天线分集增益,检测性能得到了提高。 ω=1.4 时,SOR-MCMC算法检测性能最好,后面的仿真中都取该值;ω=1 即传统MCMC算法。

图2 是收发天线数相同时, 不同天线配置下几种算法的性能比较, 采用4QAM调制,ω=1.4,T=20。由图可知, 随着天线数的增加,MCMC和SOR-MCMC算法的性能都有提高,体现了大规模MIMO的特性。高信噪比时,传统MCMC算法由于陷入问题影响了检测性能, 而SOR - MCMC算法克服了这个问题。在16 × 16 / 32 × 32 / 64 ×64 等几种天线配置下, SOR - MCMC的性能都优于传统MCMC方法, 且随着信噪比的增加, 不断逼近单用户的MFD算法。

大规模MIMO的实际应用中, 基站天线数一般远大于用户数, 下面考虑收发天线数不相等时的算法性能。基站天线数N=32、用户数K分别为8/16/32 时各算法性能如图3,仍采用4QAM调制,ω=1.4,T=20。在图3 所示的几种天线配置下,SOR -MCMC的性能都优于传统MCMC方法, 还解决了传统MCMC方法高信噪比时的陷入问题。 β 较大时,SOR-MCMC算法性能随着SNR的增大不断逼近MFD。由图3 可知,在天线配置为8×32 时,SOR - MCMC算法误码率达到10-3所需的信噪比与MFD相比只相差0.4 d B。

采用16QAM调制, 其余参数不变,SOR-MCMC算法性能仿真结果如图4。可以看到, 随着信噪比的增加,SOR-MCMC算法检测性能逼近MFD,即该算法在高阶调制下仍然有效。

4 结论

MCMC 篇4

在地质流体力学中,以孔隙度和渗透率作为参数的流量概率密度分布,往往呈现出较为复杂的函数结构,很难用常用的经典分布(如标准正态分布,均匀分布等)予以描述。因此,需要从孔隙度、渗透率及油井产量所构成的复杂参数/分布结构中获得必要的模拟样本。该文将马尔可夫链蒙特卡罗方法(Markov chain Monte Carlo method)应用于多孔介质流体预测中,介绍了MCMC方法的由来及其基本概念,阐述并依据实例展示了MCMC及其相关算法在石油产量预测中的运用,同时也对如何进一步提高预测成效提出了一定设想。

1 马尔可夫链蒙特卡罗方法

1.1 马尔可夫链(Markov Chain)

马尔可夫链[1]指的是时间和状态都离散的马尔可夫过程,简记为它是随机变量的一个数列。这些变量的范围,即他们所有可能取值的集合,被称为“状态空间”,Xn的值表示在时间n的状态。如果Xn + 1对于过去状态的条件概率分布仅是Xn的一个函数,则

一个马尔可夫链模型可以表示为(S,P,Q)。其中S表示系统所有可能的状态组成的非空的状态集,可以是有限的、可列的集合或任意非空集。P表示系统的状态转移概率矩阵,Pij表示系统在时刻t处于状态i并且在下一时刻t+1处于状态j的概率。Q指的是系统的初始概率分布。

1.2 马尔可夫链蒙特卡罗方法

马尔可夫链蒙特卡罗方法(Markov chain Monte Carlo method),简称MCMC方法,是以动态构造马尔可夫链为基础,通过遍历性约束来实现模拟目标分布的一类随机模拟方法[2]。在统计学上,MCMC方法是指从已知概率分布中采样,进而构造马尔可夫链,将其均衡分布作为所需的理想分布的算法。经典的马尔可夫链蒙特卡罗采样方法有Gibbs采样法,Metropolis采样法等。

1.3 Metropolis—Hastings算法[3]

在满足上述三个条件的情况下,Q(x|y)的选取是任意的,在实际运用中,通常选择以Y为中心的正态分布,这样做的好处是,在模拟取样过程中,越接近于上一个接受样本Y的X越容易成为随机游走的新样本。算法的具体步骤如下:

首先选择任意初始值X0作为第一个样本。对于任一状态n=0,1,2,3...假设状态n时当前样本值为xn,然后进行如下步骤:

1) 从提议函数Q(x'|xn)中产生一个候选样本x';

2) 计算接受概率α=min(1,π(x')/π(xn));

3) 以概率α置状态n+1的样本值为x';以概率1α置其为xn。

如通过相关物理地质分析可得到以孔隙率和渗透率等作为参数的油井产量的概率密度函数π(x),其函数形式由于地质构造公式的制约将会相当复杂,导致无法直接获取产量的样本,通过上述算法,我们可以任意选取初始值开始取样,一段时间后(初始化阶段结束后),取样值将会依概率收敛于目标密度函数,此时可将初始化阶段结束后得到的样本作为一组从原概率密度分布得到的近似样本。在实际运用中,由于孔隙率渗透率等地质参数是未知的,需要对其进行贝叶斯估计,其过程中同样将运用MCMC方法。

2 MCMC方法应用实例

在V. Ginting等人的研究中,应用MCMC方法对石油储层产量曲线进行了预测[4]。在本案例中,首先通过部分已知的流量数据来进行渗透率和孔隙度的贝叶斯估计,并由已知的相关分布综合资料来创建地下行为模拟。再使用这些数据运行流体模拟软件得出产量曲线。具体做法如下:

石油储层的产量曲线预测包含两个步骤:描述和预测。一般我们无法获得具体的渗透率和孔隙度的资料,但是在MCMC方法的支持下我们可以通过贝叶斯估计得到。首先依据常用假设对基于渗透率和孔隙率的分布函数进行预处理,通过运用Karhunen—Loeve(简称KL)展开降低参数维度,得到包含渗透率和孔隙率的参数向量ψ,为了得到ψ的贝叶斯估计,需要从关于渗透率和孔隙率的后验概率分布中采样,即从P(ψ|Fm)中采样,这里的Fm代表了部分已知流量数据。根据贝叶斯定理,我们知道P(ψ|Fm)∝P(Fm|ψ)P(ψ),这里P(Fm|ψ)代表了给定参数条件下的流量似然分布,而P(ψ)则代表了渗透率和孔隙率的先验概率分布。这样得到的后验概率分布形式非常复杂,因此需用Metropolis—Hastings MCMC方法进行取样。在所得样本的基础上计算后验均值、方差等特征参数作为渗透率和孔隙率的贝叶斯估计进而利用此估计创建如图1所示的地下行为模拟。

图1 (1)和(2)表示地下的渗透率和孔隙度分布;(3)和(4)分别表示t=0.4PVI和t=1.4PVI时的水饱和度。注水井(左下圆点)和生产井(右上圆点corner well和右中圆点center well)的位置如图所示。MCMCMCMC

在方法预测中,至关重要的一步是通过贝叶斯方法获取渗透率和孔隙率后验概率分布样本数据,使用这些样本数据来重建地下流体模型,再根据地下模型的渗透率和孔隙度数据来预测生产曲线的走势。我们称通过预测获得的这条曲线为预测生产曲线。图2是运用MCMC方法得出的流量预测,可以看出不论是单采样MCMC还是多采样MCMC其最终结果都跟真实参考值相当接近。

然而正如上文中提到的,实现这一方法的关键点是得到关于渗透率和孔隙率的合理的贝叶斯估计。这除了要求对MCMC及其相关算法的熟练掌握外,最重要的是对先验概率分布的合理假设。在研究中,为了得到合理的先验概率分布,提出了用KarKarhunen—Loevehunen—Loeve展开降维,然而在进行模拟运算时,为了简化模型,还是基于渗透率和孔隙率相互独立的基本假设,这也就为未来的相关研究留下了很大的提升空间。

3 关于MCMC方法的改进设想

在以上关于MCMC方法的应用实例中,为了简化模型,假设渗透率和孔隙度是不相关的,但是实际上并非如此。相关文献[5,6]表明,低渗透率低孔隙度岩石的孔隙度和渗透率之间并非相互独立。所以为了得到更精确更符合实际的预测结果,在MCMC方法中,可以在渗透率和孔隙度之间加入他们之间的关系参数,通过对这类相关系数或协方差等参数的数字估计,整个模型的精度和预测效果有望大幅度提高。除了运用Karhunen—Loeve展开的降维方法以外,还有多种统计方法可以应用于此类研究,如考虑将因子分析、主成分分析等多元统计方法,或是lasso等针对高维参数的估计方法加入研究,相关的估计与预测分析将会有进一步的发展。

4结束语

MCMC方法已广泛应用于金融、地质、环境保护、电力系统等领域的研究[7-10],应用MCMC方法研究期货市场流动、地震参数、水污染溯源、大型电力系统可靠性评估等方面的工作日趋活跃,MCMC方法作为统计学的一种重要方法,将继续在其他各学科中得到应用与发展。

MCMC 篇5

SV模型在金融领域中有着广泛的用途, 因此大多数的学者从不同的角度出发, 提出多种SV模型及其相应的估计方法。本文中带跳的随机波动模型是SV模型的改进型, 能更好的拟合股票的价格波动。但是, 也正是因为SV模型中包含着潜在变量, 涉及的似然函数和无条件矩要通过高维积分来计算, 极大似然法不能直接求解。Jacquier E, Polson N G和Rossi PE.于1994年在Journal of Business&Economic Statistics期刊上发表的文章中开创了一种分析条件方差自回归的SV模型的新技术。其中用到了Metropolis算法来构建马尔科夫链模拟工具, 并对贝叶斯和最大似然估计的性能上进行了比较, 得出了基于贝叶斯估计的马尔科夫链蒙特卡罗方法 (MCMC) 在随机波动模型分析上更有效的结论。故本文运用MCMC方法对带跳的随机波动模型进行参数估计并对上证指数进行实证分析。

1带跳的SV模型

1.1模型的贝叶斯推断

应用MCMC方法对模型进行参数估计的基础贝叶斯理论, 首先通过贝叶斯理论求得各个参数和缺失变量的后验分布密度。然后对参数样本进行Gibbs抽样或MH抽样, 最终得到参数的估计值。本文中对各个参数进行了21000次迭代, 去除初始的1000次迭代, 保证各个参数的收敛性。考虑到各个参数在数值上可能出现的偶然性, 本文中各个参数的估计值为各个参数20000次迭代的均值。

模型的联合分布密度函数为:

模型各个参数的先验分布假定为:µ~N (0, 5) , µh~N (0, 5) , Jt~Bern (λ) , Zt~N (uj, σj) , φh~N (0.95, 1) , σh2~IG (10, 0.19) λ~B (20, 30) , σj2~IG (5, 20) , µj~N (0, 0.1) 。根据参数的先验分布及似然函数, 可得出各个参数的后验分布。

1.2参数后验分布密度函数

对于后验分布密度函数为已知标准形式可直接运用Gibbs抽样;对于后验分布密度函数为非标准形式的, 可以进行Metropolis-Hastings抽样, 选取合适的建议分布, 计算接受概率, 并抽取样本。

首先求得参数的后验分布, 设θ代表所有参数θ= (µ, µh, φh, λ, σh2, µj, σj2) ,

µ的先验分布为µ~N (0, 5) ,

Zt的后验分布, 状态变量Zt的后验分布分两种状况, 当Jt=1时, Zt~N (α1, β1) ,

对各个参数进行Gibbs抽样;参数中φh的后验分布是非标准的, 故用MH方法对φh进行抽样。

2实证结果分析

本文对2005-2014年10年的上证指数的2345条日收益率数据进行实证分析。图1为上证指数的收益率时序图。

从表1中可以看出上证指数收益呈左偏形态 (偏度<0) , 且具有尖峰厚尾特征 (峰度>3) , 可以采用SV类模型建模。

应用MCMC方法, 对带跳的SV模型进行参数估计, 得出以下结果:

通过表2中各个参数的标准差可得出运用MCMC方法估计得出的参数中µh的标准差较大, 波动幅度较大, 其他各个参数标准差都较小, 波动幅度小, 较为稳定。

由图3与图1对比可得出, Jt=1时主要分布在2006年4月以后, 此时股市开始有小幅震荡, 2007年和2008年股市的震荡幅度最大, 上证指数波动幅度也十分剧烈, 此后一直震荡不断。图2中h的估计值图像也很好的描述了y值的改变, 有着实质性的改变, 从2006年4月开始, h值逐渐升高, 到2008年1月时到达最高, 也是y值震荡幅度最大的时候, 然后逐渐下降, 之后持续长期小幅波动。上证指数的收益波动具有较强的持续性, 并随h的估计值的波动而波动, 基于MCMC方法的带跳随机波动模型能够很好的模拟上证指数的收益波动。

3结语

本文对带跳的随机波动模型进行了贝叶斯分析, 设定模型参数的先验分布, 构造了基于Gibbs抽样的MCMC数值计算过程, 并对上证指数进行实证研究。研究结果表明, 基于MCMC方法的对带跳的随机波动模型能够很好的模拟我国股市的波动, 并且证明了我国股市具有较强的波动持续性。

摘要：针对股票市场波动性表现出的时变特点与“集聚效应”, 本文对带跳的随机波动模型进行研究。应用MCMC方法对模型的参数、随机波动率及跳时间进行估计, 并对上证指数进行实证分析。

关键词：MCMC方法,带跳的随机波动模型

参考文献

MCMC 篇6

视频中运动目标的跟踪技术在军事、民用交通检测等各个方面有着越来越多的应用。许多情况下,因为人们真正关注的是运动目标的信息,比如交通流量的检测、商场里的视频监控,以及军事上的入侵检测等等,这使得运动目标的跟踪成为视频处理的热点之一。然而,这个问题也是视频处理的难点之一,尤其许多视频是在摄像机运动的情况下拍摄的。

经过国内外长期的努力,视频运动目标检测与跟踪方面的研究已经取得了一定的进展。处理跟踪问题,常用且实用的方法之一有根据贝叶斯滤波得到的卡尔曼滤波法[1,2,3],这种方法在一定程度上解决了运动方式较简单情况下的跟踪问题,但该模型对视频及噪声的要求苛刻,适用性不强,无法胜任存在目标变形、遮挡和多目标等场合,且存在精度不够的现象。粒子滤波法[4,5,6,7]是贝叶斯估计的另一种应用,有效解决了适应性问题,这种方法通过随机采样来获得粒子,并用粒子之和来逼近连续的概率密度函数值。该方法得到了比较广泛的应用,但存在粒子退化的问题,目标丢失率也较高。MCMC方法[8,9,10,11,12,13]继承了粒子滤波的优点,用MCMC抽样取代传统的重点抽样,克服了粒子退化的问题,得到最广泛的应用。尤其是AMCMC(Added Markov Chain Monte Carlo)算法[14,15],作为一种改进的MCMC算法,不仅解决了MCMC方法不能跟踪个数变化目标的问题,也简化了其他改进MCMC算法的方法和步骤。

以上方法都一定程度上实现了视频运动目标跟踪,但是在背景运动的情况下,尤其是摄像机运动剧烈时,丢失目标的概率仍然很高。在本文中,通过结合代表摄像头运动的FOE[16,17,18,19],修正因摄像头运动引进的误差,有效解决了这个问题。

1MCMC粒子滤波跟踪算法

1.1贝叶斯统计算法

贝叶斯方法将未知参数看作是随机变量,使用先验概率和当前观测信息计算后验概率。其实质是用已有的信息来构造系统状态变量的后验概率密度,即通过得到的观测序列Z1:t来计算实际状态Xt。那么我们会得到不同的置信度p(Xt|Z1:t),从而得到Xt的最优估计。根据观察向量和已知的状态向量推导现在状态的预测方程为:

p(Xt|Z1:t-1)=∫p(Xt|Xt-1)p(Xt-1|Z1:t-1)dXt-1 (1)

根据贝叶斯理论,其更新方程为:

$p (X_{t} | Ζ_{1 : t}) = \frac{p (Ζ_{t} | X_{t}) p (X_{t} | Ζ_{1 : t - 1})}{p (Ζ_{t} | Ζ_{1 : t - 1})}$ (2)

由于观察过程是相互独立的,式(2)中p(Zt|Z1:t-1)可以归一化为一个常数,我们用C来代替它,然后将式(1)带入式(2),得到:

1.2基本粒子滤波

式(3)用于跟踪时,主要有两种情况,一种是在特殊确定的条件下,比如令所有的信号过程是线性和高斯的,我们可以直接求出最优估计,显然这种情况就是我们常见的卡尔曼滤波法。然而大部分信号过程并非是高斯白噪声作用下的线性过程,显然这种方法有很大的局限性。

粒子滤波法解决了这一问题,最先由Gordon提出,它为离散时间的递归滤波问题提供了一种近似的贝叶斯解决方法,其基本思想是构造一个基于样本的后验概率密度函数。这类方法通常称为粒子滤波器或序贯蒙特卡洛方法,它通过一组带权w $_{t}^{(i)}$ 的交互采样点X $_{t}^{(i)}$ 来近似求取式(3),即:

基本粒子滤波算法的一个主要问题是粒子退化,即经过几次迭代之后,差不多所有的粒子都具有负的权值,当大多数粒子的权重较小而只有少数粒子占有绝对性的权重时,经重采样后得到大量完全相同的粒子,粒子的多样性消失。

1.3MCMC算法

MCMC是在普通粒子滤波的基础上建立起来的,也是来自贝叶斯估计方程。对式(3),用马尔科夫链离散抽样来近似递归的贝叶斯滤波器的概率分布,可以得到:

式(3)变为:

MCMC算法是在目标数量固定不变的情况下,在贝叶斯多目标跟踪粒子滤波器(式(4))中,用MCMC抽样代替传统的重点抽样,产生N个不加权的粒子代替时刻的后验概率,即:

常用的MCMC实现方法是Metropolis-Hastings算法,其具体步骤为:

(1) 初始化马尔可夫链Xt(0),并对其迭代抽样,生成Nm+Nb个抽样,其中Nb是冗余的抽样个数,Nm是马尔可夫链达到收敛所需要的抽样个数。

(2) 从基于抽样n-1的提议分布q(X $_{t}^{*}$ /X $_{t}^{(n - 1)}$ )中产生一个提议抽样X $_{t}^{*}$ 。

(3) 计算接收率:

(4) 接收/拒绝判断。如果α≥1,接受提议抽样X $_{t}^{*}$ ,并把它加进马尔可夫链,即X $_{t}^{(n)}$ =X $_{t}^{*}$ ,否则提议抽样以概率α被接收。如果提议被拒绝,把前一个抽样作为当前抽样,即X $_{t}^{(n)}$ =X $_{t}^{(n - 1)}$ 。

2相机聚焦点的计算和作用

FOE 是指摄像机仅仅平移运动时,作为背景的像素点光流的交叉点。它是由于摄像机的运动造成的,代表着摄像机平移运动的方向,除了用来表示摄像头的运动方向外,它可以用来推算摄像头的运动方式,背景的运动方式等。

FOE的计算主要是利用光流场来完成的。由于密集光流场的计算复杂,并且在纹理稀疏的区域准确率明显不足,所以一般利用稀疏光流场来估算FOE。在取得角点和交点后,我们利用Lucas-Kanade算法来解决光流的计算问题。

在获得光流场后,以下推导如何计算FOE。在图1(a)的三维坐标中,Z轴代表着相机的光轴方向,也就是与成像平面垂直:X_Y轴与像平面是平行的。点P(X,Y,Z)表示三维空间中的一个点, p(x, y)是它在图像上的投影点。如果摄像头只存在平移运动,那么p点在图像平面上的运动场(光流场)可以表示为:

${\begin{cases} v x = \frac{Τ_{z x} - Τ_{x} f}{Ζ} \\ v y = \frac{Τ_{z y} - Τ_{y} f}{Ζ} \end{cases}$

(9)

式中,Tx、Ty和Tz分别表示三维平移运动在三个方向的分量,而f是摄像机的焦距。一般FOE定义为:

${\begin{cases} x_{0} = \frac{f Τ_{x}}{Τ_{z}} \\ y_{0} = \frac{f Τ_{y}}{Τ_{z}} \end{cases}$

(10)

于是式(9)可以表示成:

${\begin{cases} v_{x} = (x - x_{0}) \frac{Τ_{z}}{Ζ} \\ v_{y} = (y - y_{0}) \frac{Τ_{z}}{Ζ} \end{cases}$

(11)

这样还可以用光流分量来表示FOE坐标,同时消去Z:

vyx0-vxy0+vxy-vyx=0 (12)

写成矩阵形式:

$(\begin{matrix} v y_{1} & - v x_{1} \\ v y_{2} & - v x_{2} \\ ⋮ & ⋮ \\ v y_{Ν} & - v x_{Ν} \end{matrix}) (\begin{array}{l} x_{0} \\ y_{0} \end{array}) = (\begin{matrix} v y_{1} x_{1} & - v x_{1} y_{1} \\ v y_{2} x_{2} & - v x_{2} y_{2} \\ ⋮ & ⋮ \\ v y_{Ν} x_{Ν} & - v x_{Ν} y_{Ν} \end{matrix})$

(13)

通过计算稀疏光流场,最后作奇异值分解(SVD),就可求得FOE坐标。在实际应用中并非所有求得的光流场都能用来计算FOE,因为只有背景的部分才能代表摄像头的运动,所以我们只有选择合适的门限值,去除太大和太小的光流值,才能获得较准确的FOE坐标。图2给出一个求得FOE的示例,其中中心处大的圆点代表FOE位置。

3采用FOE的MCMC跟踪算法

我们知道,直接利用MCMC算法跟踪视频中运动目标有不错的效果,但是被用来跟踪基于运动摄像头拍摄的视频,往往会非常频繁出现目标丢失或者抖动混乱的情形。究其原因,是由于摄像头运动的加入,目标前后两帧中的位置偏离太大,超出了现有MCMC的搜寻范围,对抽样次数的需求量势必增多。如果该要求不能满足,经过几次误差积累,将导致目标丢失现象的发生。当然我们可以通过增加抽样次数,加长马尔科夫链来改善跟踪效果。但MCMC算法的计算量,主要集中在抽样的过程中,所以这样会使计算量大大增加。尤其是摄像头的运动速度是不定的,我们无法估算应该增加多少个抽样,才能正确地逼近目标位置的分布。若摄像头停止,增加的抽样次数被浪费掉,徒劳增加了抽样次数计算量。因为抽样次数在整个过程中一般是固定的,所以通过增加冗余抽样的方法有其局限性。至于抖动的产生,是因为目标在视频中的消失是需要一个过程的,当目标离开场景时,必须有一个判断的过程。所以目标离开,对其进行的采样不会立刻停止,这就使得目标采样会在目标原来位置进行一段时间的抖动。尤其是当两个目标非常近的时候,无论采用怎样精确的算法,误跟踪是不可完全避免的。

针对上述情况,当摄像头做平移运动时,我们首先计算背景角点和交叉点处的稀疏光流场,然后计算FOE。利用FOE的坐标位置,可以对目标在下一帧中的基本位置进行预测。并且,当目标位置速度相对于摄像头速度较小时,甚至可以忽略目标的运动,把目标的运动近似看成全部由摄像头的运动造成的。这样就可以对目标的位置有个基本的判断。

图3表示了预测的基本模型,假设目标在与摄像头运动方向垂直的平面上与光轴的距离为x,与摄像头的垂直距离为d,其在成像面上距离FOE的距离为x1。下一帧中,摄像头向后运动Δd,在极短时间内,假设光轴方向是不变化的,这样目标近似向前运动了Δd,在成像面上到FOE的距离为x2。由成像原理,可近似得到:

$\frac{x}{x_{1}} = \frac{d}{f}$ (14)

$\frac{x}{x_{2}} = \frac{d + Δ d}{f}$ (15)

假设摄像头即时速度为v,一帧的间隔时间为Δt,焦距为f。经过变换可得:

$Δ x = (\frac{v Δ t}{d + v Δ t}) x_{1}$ (16)

同理,若是摄像头向前运动,有:

$Δ x = (\frac{v Δ t}{d - v Δ t}) x_{1}$ (17)

以上情况是在假设目标不运动的情况下的运动模型,大部分情况是单摄像头且目标运动的,一般无法判断距离d,摄像头的速度v也是不易精确获得。本文中采用比较直观的方式:设定一个阈值,在目标物面积小于阈值的时候,这代表目标距离摄像头太远,其造成的影响可忽略;当面积大于阈值时,再利用目标位置与FOE的位置的距离和目标速度的估算,对目标位置进行修正。这样,二维坐标预测方程可简单表述成:

${\begin{cases} X_{t} = X_{t - 1} + a (X_{t - 1} - X_{f o e}) + Δ X \\ Y_{t} = Y_{t - 1} + a (Y_{t - 1} - Y_{f o e}) + Δ Y \end{cases}$

(18)

式中a代表着摄像头与目标的距离以及载体运动速度对位置的综合影响,在这是一个经验值。其正负表示了摄像头是向前还是向后运动。如果是前进中拍摄的,a为正值,代表摄像头运动促使目标发散;若摄像头向后运动,a为负值,代表摄像头运动促使目标像图像中心汇聚。整个式子表示由摄像头运动所引起的目标运动,处于和摄像头同一直线的目标所受影响较小,而在远离坐标点,也就是和摄像头所在直线较远的位置,受到的影响较大。通过这样一个修正,能对目标的基本运动分量进行预测,从而有效修正了目标位置。并且,预测值的加入,基本规定了运动目标的消失方向,有效减少了目标消失后,四处搜寻目标的现象。

Metropolis-Hastings算法,是MCMC实现的常用方法之一,该方法和FOE的结合方法,如图4所示。其右面代表了普通Metropolis-Hastings算法,左边代表FOE的计算过程。FOE与MCMC的结合,主要改变在抽样这一过程。抽样时,在预测位置处进行抽样。实际操作中,因为预测是不够精确的,所以加入一个随机函数,在一次抽样中,只对一部分的目标进行预测。以使结果更加准确,避免因为不精确预测引进误差。先根据摄像头的运动(FOE),对目标的位置进行一个简单的线性预测,然后在此基础上进行抽样。这样就消除了一部分因为摄像头运动带来的误差,有效地减少了目标跟丢的现象。

4试验结果

实验视频模仿海底潜艇拍摄海底生物,试验环境是PC机(Pentium 2.5GHz,内存1.99G)利用VS2005开发。视频图像的分辨率是320×240,是水下潜艇在高速运行的情况下拍摄的海底鱼类、海蜇以及其他生物,整个过程中潜艇向前运动。

图5和图6分别显示了未改进MCMC和改进后MCMC方法在同样条件下跟踪三个目标的过程。

由结果可见,在跟踪相同目标的情况下,未改进的方法出现了跟丢的情况,这是因为该视频是在摄像头前进的情况下拍摄的,目标运动幅度较大,可能会超过MCMC的搜寻范围。而与单纯的MCMC算法相比,本算法对摄像头的向前运动进行了预算,通过FOE估计目标在场景将要出现的位置,并在此位置处进行搜寻,得到了更为精确的跟踪结果。

在图5中,当有目标离开场景时,其并未引进标定框四处搜寻的现象。这说明本算法在一定程度上消减了目标离开场景时的抖动。除了上述效果,我们还发现在实现相同跟踪目的的情况下,能以较少的抽样次数,实现相同的目的,使程序运算速度加快。

5结语

本文介绍了一种在摄像头运动的情况下视频中运动目标跟踪的改进方法,通过计算代表摄像头运动的FOE和目标的平移分量速度,对运动目标的基本运动分量进行估计,达到对运动目标位置的修正与预测。在与AMCMC结合使用时,减少了目标丢失的数目,并有效消除了目标在边缘处抖动的情况。假如在能获得目标到摄像头距离的情况下,本算法能更精确地判定目标位置,进一步提高跟踪的准确性及稳定性。

MCMC 篇7

粒子滤波器是一种基于蒙特卡罗和递推贝叶斯估计的算法。基本思想是在状态空间中用一系列加权样本来近似后验概率密度分布,以样本均值代替积分运算。原则上,粒子滤波器可以实现任意状态的估计,尤其在卡尔曼滤波器和扩展卡尔曼滤波器失效的非线性或非高斯状态估计中表现出优越的性能。粒子滤波器己经在目标跟踪、故障诊断、导航定位、无线通信和信息管理[1]等领域得到广泛应用。

1993年,由Gordon等提出的一种新的基于SIS的非线性滤波方法[2],奠定了粒子滤波算法的基础。退化现象是粒子滤波普遍存在的一个问题,即经过几次迭代之后,差不多所有的粒子都具有负的权值,这就意味着大量的计算都用来更新粒子,而这些粒子对逼近状态概率密度函数的贡献几乎为零。显然,退化问题在粒子滤波中是不被期望的,减小这种作用的一个简单方法就是采用非常大的样本容量,但这通常是不实际的,因为这样必然影响计算的实时性。一般粒子滤波的执行效率取决于粒子的数量,而粒子的数量由状态方程的维数,先验概率密度函数和重要密度函数的相似度以及迭代次数决定[3]。因此,解决退化问题的正确方法是引导粒子向高似然区移动。为此,可以考虑的两种办法是:优选重要性密度函数和重采样。

随着对粒子滤波算法研究的深入,不断地有学者对算法的不足之处进行完善,并将其应用到越来越广泛的领域。本文首先研究了粒子滤波的改进算法:U粒子滤波算法和MCMC粒子滤波算法。针对粒子滤波算法存在退化现象的缺点,本文将两种算法结合在一起,以提高算法的性能。为了改善算法运算量大的缺点,对UPFMCMC算法进行了进一步的改进,以提高运行效率,仿真实验验证了这种改进方法的优越性。

1 U粒子滤波

在许多实际应用问题中,当状态方程或测量方程为非线性并且噪声为非高斯情况时,滤波问题也表现为非线性。由于近似非线性函数的概率分布比近似非线性函数更容易,使用采样方法近似非线性分布来解决问题的途径得到了人们的广泛关注。文献[4]提出了采用U卡尔曼滤波(UKF)[5]方法对非线性问题进行滤波估计。UKF是一类采用策略逼近非线性分布的方法。UKF以Unscented变换(简称UT变换)为基础,采用卡尔曼滤波(Kalman Filter,KF)框架,具体采样形式为确定性采样,而非粒子滤波的随机采样,从而避免了粒子滤波的粒子退化问题。

UKF能够很容易地将最新观测量引入状态估计,但是,它要求对状态分布进行高斯假设。粒子滤波虽然适用于任何分布情况的状态估计问题,但基本粒子滤波算法很难把最新观测量引入重要分布函数。为了更好地利用UKF和粒子滤波算法各自的优点,取长补短,将UKF算法和粒子滤波算法相结合,形成U粒子滤波算法(UPF)[6]。

UPF算法的核心思想是在粒子滤波的基础上利用UKF得到比普通PF更好的重要性分布函数,即:

$π (x_{k}^{(i)} | x_{0 : k - 1}^{(i)}, z_{0 : k}) = Ν (\bar{x}_{k}^{(i)}, Ρ_{k}^{(i)})$

式中 $\bar{x}_{k}^{(i)}$ 和P $_{k}^{(i)}$ 是通过UKF算法得到的均值和协方差。利用UKF算法产生UPF算法所需要的重要性密度函数,二者有效结合便是UPF算法。

UPF算法的具体实现步骤如下:

Step 1:初始采样,x $_{0}^{(i)}$ ～p(x0),设置w $_{0}^{(i)}$ =1/N,i=1,2,…,N;

Step 2:对时刻k(k∈N),使用UKF算法求粒子集{x $_{k}^{(i)}$ } $_{i = 1}^{Ν}$ 的均值 $\bar{x}_{k}^{(i)} = E (x_{k}^{(i)})$ ,方差P $_{k}^{(i)}$ ;

Step 3:从重要性密度函数, $π (x_{k}^{(i)} | x_{0 : k - 1}^{(i)}, z_{0 : k}) = Ν (\bar{x}_{k}^{(i)}, Ρ_{k}^{(i)})$ 中采样得到粒子x $_{k}^{(i)}$ ,i=1,2,…,N;

Step 4:根据式(1)得到粒子权值w $_{k}^{* (i)}$ :

$w_{k}^{* (i)} = w_{k - 1}^{* (i)} \frac{p (y_{k} | x_{k}^{(i)}) p (x_{k}^{(i)} | x_{k}^{(i)})}{π (x_{k}^{(i)} | x_{0 : k - 1}^{(i)}, z_{0 : k})} (1)$

Step 5:根据式(2)计算,得到归一化权值w $_{k}^{(i)}$ :

$w_{k}^{(i)} = w_{k}^{* (i)} / \sum_{i = 1}^{Ν} w_{k}^{* (i)} (2)$

Step 6:根据式(3)更新状态估计值:

$\hat{x}_{k} = \sum_{i = 1}^{Ν} x_{k}^{(i)} w_{k}^{(i)} (3)$

2 MCMC粒子滤波

粒子滤波算法中的重采样[7]步骤在一定程度上可以减小退化问题,但因此又会引来粒子耗尽问题(采样贫乏),即经过几次迭代后,粒子中将包含许多重复点,这在估计过程中会由于粒子失去多样性而使后验概率的离散逼近变得不精确,从而影响粒子滤波器的估计性能。在动态噪声较小时,这一问题尤为严重。为了解决粒子耗尽问题所带来的副作用,恢复粒子的多样性,在重采样步骤后引入马尔可夫链蒙特卡罗移动步骤(MCMC)[8],使得粒子分布更加合理。

MCMC的基本思想[9]是:如果粒子的分布服从重要性函数 $q (\hat{x}_{0 : t} | z_{1 : t})$ ,则可以应用一个马尔科夫链转换核 $κ (x_{0 : t} | \hat{x}_{0 : t})$ 和后验概率密度函数p(x0:t|z1:t),使得:

$\int \begin{matrix} κ ( \end{matrix} x_{0 : t} | \hat{x}_{0 : t}) q (\hat{x}_{0 : t} | z_{1 : t}) d \hat{x}_{0 : t} = p (x_{0 : t} | z_{1 : t})$

由此产生的新粒子仍服从状态的后验概率分布,但新的粒子群更加逼近真实目标分布。MCMC转移步骤实际就是从有限混合分布 $\frac{1}{Ν} \sum_{i = 1}^{Ν} κ (x_{0 : t} - \hat{x}_{0 : t}^{(i)})$ 中进行采样来增加粒子的多样性。

MCMC转移有很多实现算法(Gibbs采样,Metropolis Hastings算法等)。本文采用MH算法来实现MCMC转移,其步骤如下:

Step 1:产生一个在[0,1]区间上服从均匀分布的随机数v;

Step 2:从重要性函数 $q (\overset{⌒}{x_{k}} | x_{0 : k - 1}^{(i)}, z_{1 : k})$ 中采样产生备用粒子x $_{k}^{* (i)}$ ,i=1,2,…,N;

Step 3:如果v≤min{1,α},则接受x $_{k}^{* (i)}$ ,这时 $x_{0 : k}^{(i)} = (\overset{⌢}{x}_{0 : k - 1}^{(i)}, x_{k}^{* (i)})$ ;否则拒绝 $x_{k}^{* (i)} ‚ x_{0 : k}^{(i)} = \overset{⌢}{x}_{0 : k}^{(i)}$ 。

其中α称为接受比率:

$α = \frac{p (z_{k} | x_{k}^{* (i)}) p (x_{k}^{*} | \overset{⌢}{x}_{k - 1}^{(i)}) q (\overset{⌢}{x}_{k} | \overset{⌢}{x}_{0 : k - 1}^{(i)}, z_{1 : k})}{p (z_{k} | \overset{⌢}{x}_{k}^{(i)}) p (\overset{⌢}{x}_{k} | \overset{⌢}{x}_{k - 1}) q (x_{k}^{* (i)} | \overset{⌢}{x}_{0 : k - 1}^{(i)}, z_{1 : k})}$

将上述算法加在粒子滤波的重采样步骤后,即是MCMC粒子滤波算法。

MCMC粒子滤波算法的基本步骤归纳如下:

Step 1:初始化,产生粒子集{x $_{0}^{(i)}$ } $_{i = 1}^{Ν}$ ,所有粒子的权值均为1/N;

Step 2:重要性采样

(1) 从重要性函数π(xk|x $_{k - 1}^{(i)}$ ,yk)中采样得到粒子集{x $_{k}^{(i)}$ } $_{i = 1}^{Ν}$ ;

(2) 根据式(1)计算k时刻的权值;

(3) 根据式(2)归一化权值;

Step 3:重采样

(1) 根据式 $Ν_{e f f} = \frac{1}{\sum_{i = 1}^{Ν} (w_{k}^{i})^{2}}$ 计算有效抽样尺度的近似值;

(2) 当 $\overset{⌢}{Ν}_{e f f}$ 小于设定的门限Nth时,进行重采样;

(3) 得到新的粒子集{x $_{0 : k}^{* (i)}$ } $_{i = 1}^{Ν}$ ,其归一化权值均为w $_{k}^{(i)}$ =1/N;

Step 4:MCMC移动

(1) 产生一个在[0,1]区间上服从均匀分布的随机数v;

(2) 从重要性函数 $q (\overset{⌢}{x}_{k} | x_{0 : k - 1}^{(i)}, z_{1 : k})$ 中采样产生备用粒子x $_{k}^{* (i)}$ ,i=1,2,…,N;

(3) 如果v≤min{1,α},则接受x $_{k}^{* (i)}$ ,这时 $x_{0 : k}^{(i)} = (\overset{⌢}{x}_{0 : k - 1}^{(i)}, x_{k}^{* (i)})$ ;否则拒绝 $x_{k}^{* (i)} ‚ x_{0 : k}^{(i)} = \overset{⌢}{x}_{0 : k}^{(i)}$ ;

Step 5:根据式(3)更新状态估计值。

3 改进的算法

对于粒子滤波退化问题的改善主要在优选重要性密度函数和重采样两方面来进行。本文将U粒子滤波算法和MCMC粒子滤波算法相结合,产生具有MCMC移动步骤的U粒子滤波(UPFMCMC)算法。在该算法中选择了UKF将最新观测量引入状态估计,以获得更好的重要性分布;重采样算法用来抑制粒子滤波的退化现象,同时在重采样步骤后引入了统计理论中的MCMC方法,增加了粒子的多样性,改善了算法的性能。

在带有MCMC移动步骤的U粒子滤波算法中,首先UPF根据前一时刻粒子及其方差确定一组Sigma点。此点集的位置和权值是由粒子的期望和方差惟一确定的,能较准确地抓住粒子概率分布的特征,然后将这些点分别代入状态方程,得到一新点集,用这些点集的加权和作为期望,用其方差的加权和作为方差,然后再用测量方程对已经求得的期望和方差进行修正,并用修正后的值作为高斯分布的期望和方差,产生一个当前时刻的粒子,因为其充分考虑当前测量值对后验概率分布的影响,提高了粒子的利用效率,但每个粒子产生阶段的计算量较大。因此本文对UPF算法做如下改进:对前一状态的估计值进行UKF,由重要性密度函数 $π (x_{k}^{(i)} | x_{0 : k - 1}^{(i)}, z_{0 : k}) = Ν (\bar{x}_{k}^{(i)}, Ρ_{k}^{(i)})$ 产生20%粒子,其余80%的粒子由先验概率分布产生。经过这样改进的算法,既考虑当前量测、先验概率对后验概率分布的影响,又减少计算量。

其次基于MCMC的粒子滤波改进算法,虽然有效地改善了粒子枯竭问题,改善了粒子滤波算法的性能,但是,随之而来又产生了计算量大大增加的问题,并且,随着粒子滤波算法中粒子数目的增加,其计算量亦会极大地增加,使得算法的实时性变得很差,甚至会导致算法的发散。在粒子滤波算法的估计过程中,引入在基本的粒子滤波算法中己经使用过的有效性尺度[10]这一概念,通常用有效性尺度来衡量粒子滤波算法的退化程度,有效性尺度越小,表示退化越严重,即多数粒子只具有很小的权值。对有效性尺度的定义如下:

$Ν_{e f f} = r o u n d (1 / \sum_{i = 1}^{Ν} (\tilde{w}_{k}^{i})^{2}) (4)$

式中round(·)表示向最近的整数取整的运算。

当有效性尺度Neff较大时,表明粒子滤波的退化现象并不严重,那些权值很小的粒子对结果的影响微乎其微,但是对它们的计算却要占用大量的计算资源,这样势必会影响粒子滤波的计算效率。所以这时去除权值很小的粒子,这样就可以不用对这些去除的粒子进行MCMC重采样步骤,大大节约计算资源,提高粒子滤波算法的计算效率;相反,当有效性尺度Neff较小时,表示退化现象比较严重,这时就需要进行重采样。本文将有效性尺度的概念引入到基于MCMC的粒子滤波算法中,有效地提高了该算法的计算效率。

综上所述,对粒子滤波的改进算法的实现步骤如下:

Step 1:初始采样,x $_{0}^{(i)}$ ～p(x0),设置w $_{0}^{(i)}$ =1/N,i=1,2,…,N;

Step 2:对时刻k(k∈N),使用UKF算法求粒子集{x $_{k}^{(i)}$ } $_{i = 1}^{Ν}$ 的均值 $\bar{x}_{k}^{(i)} = E (x_{k}^{(i)})$ ,方差P $_{k}^{(i)}$ ;

Step 3:从重要性密度函数 $π (x_{k}^{(i)} | x_{0 : k - 1}^{(i)}, z_{0 : k}) = Ν (\bar{x}_{k}^{(i)}, Ρ_{k}^{(i)})$ 中采样得到20%的粒子x $_{k}^{(i)}$ ,i=1,2,…,0.2*N;利用PF产生剩余80%的粒子;

Step 4:根据式(1)分别得到粒子权值w $_{k}^{* (i)}$ ;

Step 5:根据式(2)计算,得到归一化权值w $_{k}^{(i)}$ ;

Step 6: 根据式(4)计算有效性尺度Neff。如果Neff小于门限值Nth,则进行步骤7,否则进行步骤8;

Step 7:MCMC重采样:

(1) 产生一个在[0,1]区间上服从均匀分布的随机数v;

(2) 从重要性函数 $q (\overset{⌢}{x}_{k} | x_{0 : k - 1}^{(i)}, z_{1 : k})$ 中采样产生备用粒子x $_{k}^{* (i)}$ ,i=1,2,…,N;

其中α称为接受比率:

Step 8:根据式(3)输出状态估计值。

对于阈值的选取,根据经验,一般情况选取为Nth=N/3。

4 仿真和结果分析

本次仿真硬件环境是Intel( R) Core2 Duo CPU T7250 2. 00 GHz,2G RAM,Windows XP操作系统,软件利用Matlab R2007b编写。

仿真选取典型的单变量非平稳模型作为估计目标,来比较U粒子滤波、MCMC粒子滤波和改进的U-MCMCPF三种算法下的状态估计性能,其中时间步长取50,独立实验次数为500,粒子数目分别取500和1 000,重采样阈值Nth选为N/3。状态估计性能以500次独立实验三种算法下的粒子滤波器平均运行时间和RMSE均值作为衡量标准。

三种不同算法下的状态估计如图1所示。

三种不同算法下的粒子滤波器运行时间和RMSE均值如表1所示。

从仿真图和仿真数据中不难看出:

(1) U粒子滤波、MCMC粒子滤波和改进的U-MCMCPF三种算法下的状态估计效果很接近。500次独立实验三种算法下的RMSE均值差别很小。

(2) 从运行时间看,改进的U-MCMCPF算法要明显优于U粒子滤波、MCMC粒子滤波两种算法,相对于U粒子滤波、MCMC粒子滤波,粒子数目为500时,改进的U-MCMCPF算法下的粒子滤波器的平均运行时间分别缩短了32.7%和54.4%,在粒子数目为1 000时,则分别缩短了32.8%和46.3%。

5 结语

本文通过对MCMC移动步骤的U粒子滤波算法进行性能上的分析,并且从样本产生和对MCMC移动步骤的使用两个方面对算法进行改进,减小了算法的计算量,提高算法的运算速度。通过仿真验证了算法的有效性和实时性。

参考文献

[1]程水英,张敛云.粒子滤波评述[J].宇航学报,2008,29(4):1099-1111.

[2]GORDON N,SALMOND D.Novel approach to nonlinearand non-Gaussian Bayesian state estimation[J].Proc.ofInstitute Electric Engineering,1993,140(2):107-113.

[3]杨小君,潘良,杨彦.基于粒子滤波和检测信息的多传感器融合跟踪[J].信息与控制,2005,34(3):356-359.

[4]MERWE V D,FREITAS N D.The unscented particle filter[C]//14th Annual Neural Information Processing Systems(NIPS)Conference.Seatrle,WA,USA:IEEE,Denver,Colorado,USA:[s.l.],2000:584-590.

[5]JULIER S J,UHLMANN J K,DURRANT-WHYTEN HF.A new approach for filter nonlinear system[C]//Proc.of the 1995 American Control Conference.Seattle,WA,USA:IEEE,1995,3:1628-1632.

[6]侯晓宇.改进的粒子滤波算法在OFDM系统盲均衡中的应用研究[D].北京:北京交通大学,2007.

[7]冯弛,王萌,汲清波.粒子滤波器的重采样算法的分析与比较[J].系统仿真学报,2009,21(4):1101-1105.

[8]朱志宇.粒子滤波算法及其应用[M].北京:科学出版社,2010.

[9]SPALL J C.Estimation via Markov chain Monte Carlo[J].IEEE Control System,2003,23(2):34-45.

MCMC 篇8

Sh ibor是央行稳步推进利率市场化改革, 提高金融机构自主定价能力, 指导货币市场产品定价, 完善货币政策传导机制而培育的中国货币市场基准利率体系。自2007年央行推出shibor利率后, 在央行进8年的积极培育下, 已基本确立了shibor利率在我国货币市场的基准利率地位, 并逐渐成为传导货币政策、反映市场利率变动的重要指标, 并在市场化利率形成机制中发挥重要作用。同时, 市场上发展和创新了一大批以Shibor为基准金融产品。大量债券交易、存款、贷款、贴现和理财类产品定价逐步与Shibor挂钩。目前, 金融机构市场成员交易层面已有部分金融产品定价与Shibor相挂钩:衍生品市场利率互换、远期利率协议;货币市场同业借款、同业存款、货币互换、理财产品等;债券市场以Shibor为基准浮动利率金融债券、企业债券、企业短期融资券, 以Shibor为基准的票据转贴现、票据回购等;以Shibor为定价及交易基准的金融产品越来越多, 如债券买卖、票据贴现、个人及机构理财产品, 最终, 银行存贷款利率定价也与Shibor相挂钩。

因此, 随着越来越多地的产品与Shibor挂钩, 我们应积极培育我国商业银行金融衍生产品的自我定价能力, 让Shibor真正在我国金融自由化和利率市场化改革进程中起到货币政策利率传导的主导与核心作用。以Shibor为标的的金融衍生产品定价问题成为越来越多学者的研究对象。要做到这一点, 我们首先要研究Sbibor及其本身的利率期限结构, 对其进行利率建模, 从而能够有效刻画出其利率期限结构动态变化的特征, 进而对利率的未来变动进行科学的预测, 这样不管是对以Shibor为标的的金融衍生产品定价, 还是利率市场化条件下以Shibor作为货币市场基准利率风向标的金融风险管理都起着极其重要的作用。

二、shibor利率模型的选择

由于金融市场兴起较早而且比较完善, 国外在利率方面的研究已经发展到了相当高的水平。在利率动态模型方面, 单因素利率模型主要有Merton (1973) 、Vasicek (1997) 、CIR (Cox, Ingersoll和R oss) (1985) 、和CKLS (Ch an, Kar oly i, Longstaff和Sander s) (1992) 等。其中, CKLS模型为不同的利率期限结构建立了一个共同框架利率模型, 单其利率模型只描述了利率漂移项的均值回复特性与利率扩散项的水平效应。Merton (1976) 在标准布朗运动中加入跳跃因素来对金融市场中存在的不连续动态变化进行描述。随后也有国外学者, 如Das (2002) , Chen&Scott (2002) 等人为在模型中加入随机波动率更贴近现实, 包含了更多的利率资产波动率时变特征。

回顾传统的利率模型, 我们知道它们都假设利率服从连续扩散过程。但是大量的实证研究表明, 现实情况中利率经常会出现突然的跳跃和巨幅的震荡, 与传统假设相违背。Das把Vasieck模型扩展到跳跃扩散模型, 并通过实证表明跳跃过程能捕捉连续模型无法解释的利率特征, 认为加入随机波动了更贴近现实。Jarrow等发现跳跃行为能很好地解释“波动率微笑”之谜。Ahn和Thompson认为利率的样本轨迹存在非连续性特征, 由此提出了利率期限结构的跳跃扩散模型。Johannes提出了一个关于跳跃行为的统计检验量, 并发现三个月的短期国债利率存在显著的跳跃行为。因此, 我们可以这样认为:在shibor利率模型中加入跳跃项能够更好地刻画和解释实际中的利率波动情况。

本文借鉴潘璐 (2011) 所用的CKLS-SV-Jump模型:

其中, corr[d Z1, t, d Z2, t]=0, J1, t~N (ψ1, γ12) , J2, t~N (ψ2, γ22) , q1, t, q2, t分别服从频率为λ1*和λ2*的泊松分布, 其中, Pr{q1, t=1}=λ1Δt, Pr{q1, t=0}=1-λ1Δt, Pr{q2, t=1}=λ2Δt, Pr{q1, t=0}=1-λ2Δt。该模型在扩散模型的基础上加入了跳跃项, 其实证研究证明该模型能很好地刻画利率的均值回复、水平效应、跳跃的特征。

三、MCMC参数模拟方法

(一) 马尔可夫链模拟

在介绍马尔可夫模拟之前, 先回顾马尔可夫过程。考虑一个随机过程{Xt}, 假定每个Xt都在空间Θ上取值, 如果h>t当时, 随机过程{Xt}满足P{Xh|Xs, s≤t}=P (Xh|Xs) , 则称过程{Xt}为马尔可夫过程。即马尔可夫过程满足这样的性质:给定Xt的值, Xh (h>0) 的值不依懒于Xs (s<t) 的取值。

MCMC方法构建的基础是Cliffor d-H ammer sley理论, 该理论认为基于参数的后验信息P (θ|x, y) 和状态变量的后验信息P (x|θ, y) , 即可唯一确定联合后验分布。而MCMC方法要求从联合后验分布中中抽取样本, 因此我们根据Clifford-Hammersley理论, 可以分别从参数的后验信息p (Θ|x, y) 和状态变量的后验信息p (x|Θ, y) 中抽取样本, 通常参数和状态变量其各自的后验分布都是比较容易得到的。由贝叶斯原理, 参数的后验分布可以写成参数的先验分布、似然函数和某个常数的乘积形式。MCMC方法的核心就在于将高维度的联合分布p (Θ, x|y) 分解为低维度的p (x|Θ, y) 和p (Θ, x|y) , 然后再从这两个分布中分别取样。

由于运用MCMC参数估计方法前需要对模型进行离散化, 因此, 我们接下来对欧拉离散法进行了介绍。

(二) 欧拉离散法

接下来对CKLS-SV-Jump模型进行离散化, 可以得到:

进一步可以得到下式:

最终我们整理成下式:

因此, 我们可以将rt+Δ基于rt的条件分布和θt+Δ基于θt的条件分布表示为:

其中, {εt}~iid.N (0, 1) ;{Jt}~iid.N (0, 1) 。

(三) Gibbs抽样

Gibbs抽样方法是MCMC方法中应用最广泛的一种抽样方法, 是用来获取一系列近似等于指定的多维概率分布的观察样本的方法。Gibbs抽样是单元素Metropolis-Hastings算法, 在Gibbs抽样中, 可能难以抽取, 而Metropolis方法具有很大的灵活性, 它可取为易于抽取的分布。

在足够的燃烧期后, Gibbs序列才能收敛到一个独立于初始值的平稳分布, 能否得到这个平稳分布是Gibbs抽样的关键。由于燃烧期内产生的参数值是不平稳的, 因此要放弃燃烧期内产生的样本。如果我们进行了n次抽样, 燃烧m次, 那么就将序列X (m+1) , X (m+2) , …, X (n) 作为马尔可夫链的实现值, 进行后续的参数估计和统计描述工作。

在Gibbs抽样中, 一个比较困难的事情是判断上述迭代过程是否收敛到均衡状态。不少作者都建议用的某些特征来判断。最常用的是观察遍历均值是否收敛来判断。比如在由Gibbs抽样得到的链中, 每隔一段距离计算一次参数的遍历均值, 未使用来计算平均值的变量近似独立, 通常可每隔一段取一个样本, 当这样得到的均值稳定后, 可认为Gibbs抽样收敛。

(四) Metropolis-Hastings算法

1953年, Metr opolis等人提出了一种构造转移核的方法, 该方法此后被Hastings加以推广, 形成了Metropolis-Hastings方法。具体思路如下:

MH算法的思想是构造一个以目标分布π (x) 为极限分布的Mar k ov链。任选一个不可约转移核概率p (·, ·) 以及一个函数a (·, ·) , 0<a (·, ·) ≤1, 对任一组合 (x, x') (x≠x') , 定义

则p (x, x') 形成一个转移核。

假设卡尔可夫链在时刻t处于状态x, 即X (t) =x, 则首先由q (·|x) 产生一个潜在的转移x→x', 然后根据概率a (x, x') 决定是否转移。即在潜在转移点x'找到后, 以概率a (x, x') 接受x'作为卡尔可夫链在t+1时刻的状态值, 以1-a (x, x') 拒绝转移到x', 即马尔可夫链在下一时刻的状态值仍为x。也就是说, 如果置X (t+1) =x', 我们就称之为“接受建议”;否则就称为“拒绝建议”。在实际计算中以概率a (x, x') 接受x', 可以这样实现:从U (0, 1) 产生一个随机数u, 如果u≤a (x, x') , 则置X (t+1) =x';否则置X (t+1) =x。

因此, 在有了x'后, 我们可以从[0, 1]上均匀分布抽一个随机数u, 有

这里, 分布q (·|x) 称为建议分布。由于我们的目标是使后验分布成为平稳分布, 所以在有了q (·|·) 后, 我们要选择一个a (·, ·) 使相应的p (x, x') 以π (x) 为其平稳分布, 我们通常选择

因此也有

该式产生的马尔可夫链是可逆的, 即π (x) p (x, x') =π (x') p (x', x) , 且π (x') 是该式确定的马尔可夫链的平稳分布。

(五) MH随机游走采样法

如果M-H算法中的建议函数q (x, x') 满足1、对称性q (x, x') =q (x', x) ;2、q (x, x') 仅与x'-x有关, 即q (x, x') =q (|x'-x|) , 那么算法就演变为通常意义下的随机游动采样法。最常见的一种随机游动采样法以正态分布, 即q (x'-x) =φ (x'-x) 为建议函数。

这里是任意一种多维正态分布的密度函数, 也就是说x'~N (x, Σ) 。其中, x为当前状态值, Σ是任意的正定矩阵。在实际计算中, 由于要求Σ是正定矩阵, 所以常取为Σ=σI, σ是一个参数。这时候, 如何选取建议函数的方差σ是影响算法效率的主要因素。对于一维情况, 就是要选择一个合适的σ值, 而对于多维情况, 我们需要确定给一个合适的正定方差矩阵Σ。一般认为在模拟多维正态随机数时, 如果调整σ的值使得在整个模拟过程中建议被接受的比例为20%左右, 则将得到比较好的模拟效果。

随机游动采样法是M-H算法中较为常用的一种形式, 们已经看到对于随机游动采样法而言, 如何选择建议函数的方差成为影响算法效率的主要因素, 对于一维情况, 就是要选择一个合适的σ值;而对于d>1的多维情况, 我们需要一个合适的正定方差矩阵Σ。在选择σ的过程中, 根据关于采用多维正态作为提议函数的随机游动采样法的一些理论结果, 一般认为如果目标分布是d维的且各分量之间独立, 则调整接受概率到0.234左右将最大限度地优化随机游动采样法的效率, 并且接受概率0.234与目标分布的具体形式无关。

基于这样的想法, 我们尝试构造一种自适应算法来寻找合适的提议函数的方差。大致的思想是让接受比率落入一个可以接受的范围内, 比如区间[a-ε, a+ε]。以一维的情况为例, 先给定一个初值σ0, 用这个值产生一条链, 来估计接受比率P0。若P0在给定的区间[a-ε, a+ε]中, 则符合要求, 算法停止;反之, 则以为中心做一次随机游动, 得到一个σ1, 再计算P1是否满足要求, 以此类推, 直至找到合适的σ值。具体实现步骤如下:

1) 给定σ一个初值σ0, 置n=0;

2) 以N (Xi-1, σn) 为建议函数, 用随机游动采样法产生一条长度为m的马尔可夫链;

3) 计算第2步中产生的马尔可夫链接受新状态的比率Pn;

5) 如果n>2且|Pn-a|<|Pn-1-a|, 则置σ=σn-1;

6) 从N (σn, σ) 中产生一个随机数, 记为σn+1, 置n=n+1, 转入第2步。

其中, 步骤5的判断是为了保证接受比率Pn要优于Pn-1, 即σn至少不会比σn-1坏;第6步中σ的选取没有固定的规律可循, 主要是看对初始σ0的估计是否接近于合适的σ值。若估计σ0与σ较为接近, 则σ可选择得小一些;若估计σ0与σ相差较远, 则σ可选择得大一些;对于d>1的多维情况, 如果要确定一个任意的正定矩阵∑, 则要确定d2个参量, 显得有些麻烦。更多情形下, 我们取∑=σI。这样一来要确定的参数就只有一个, 同时上述算法也可以用于多维情形, 只需在6步中“以N (Xi-1, σn) 为建议函数”改为“以N (Xi-1, fnJ) 为建议函数”即可。

四、模型参数估计结果

(一) 模型的先验分布

贝叶斯学派认为, 在进行观察以获得样本之前, 人们对θ也会有一些知识。因为是在试验观察之前, 故称之为先验知识。因此, 贝叶斯派认为, 应该把θ看作是随机变量。θ的分布函数记为H (θ) , θ的密度函数记为h (θ) , 分别称为先验分布函数和先验密度函数, 两者合称为先验分布。在MCMC方法中, 待估参数的先验分布的设定将会对该方法运行效率产生很大影响。本文为了方便起见, 参照了Siddhartha Chib, Federico Nardari和Neil Shephard (2002) 和潘璐 (2011) 设定的参数先验分布。根据共轭理论, 参数有相同的后验分布和先验分布。

(二) 抽样方法

波动率状态变量的完全联合后验分布可以表示为:

其中,

所以, 我们得到:

因为上述分布是不可识别的, 所以我们无法使用Gibbs抽样方法, 而需要使用Metropolis抽样方法进行抽样。而且由于它的复杂形式, 在没有采取近似处理的情况下, 用但不Metropolis方法是无法直接进行抽样的。因此, 根据Jacquier, Polson和Rossi (1994) 的研究成果, 我们采用随机游走Metropolis算法。因为我们可以对满条件分布p (r|rt-Δ, rt+Δ, Θ, r) , 尤其是它的尾部进行很好的近似。在上式中, 第一项是对数正态分布, 第二项是倒Gamma分布, 所以满条件分布p (r|rt-Δ, rt+Δ, Θ, r) 的建议分布可以很好地由某个合适的倒Gamma分布来近似。我们设这个建议分布为q (θt) , 真正的条件密度函数为π (rt) =p (r|rt-Δ, rt+Δ, Θ, r) 。

具体地, 我们将CKLS-SV-Jump模型的MCMC算法分解成以下步骤:

……

第十二步:利用Metropolis-Hastings自适应算法从p (r|rt-Δ, rt+Δ, Θ, r) 中抽取θt。

其中第十二步又可以分为以下几步:

(1) 给θ设定一个初始值θ0, 置t=0;

(2) 以N (rt-Δ, θt) 为建议分布, 用随机移动M-H算法产生一条长度为的马尔可夫链;

到此为止, 完成第一轮迭代过程。接下来, 利用新的参数作为初始值, 重复以上步骤, 得到更新参数, 重复前面的步骤次, 直至各参数和波动率潜在变量θt的模拟轨迹收敛为止, 得到一系列的随机抽取。我们可以认为 (αm*, βm*, ρm, ψ1, m, ψ2, m, γ21m, γ22m, λ1, m, λ2, m, υm, φm, τm2, θm) 渐进等价于参数估计值, 然后把收敛时的模拟值作为对各参数和潜在变量θt的估计值。

(三) 数据选择

【MCMC】推荐阅读：

>> 查看更多相关文档

MCMC

MCMC 篇1

3.1 仿真1算法的收敛性仿真

3.2 仿真2算法分类性能比较

3.3 仿真3运算量比较

MCMC 篇2

MCMC 篇3

MCMC 篇4

MCMC 篇5

MCMC 篇6

MCMC 篇7

MCMC 篇8

本站热搜

相关推荐