MCMC方法

2024-10-31

MCMC方法（精选8篇）

MCMC方法篇1

流体在多孔介质中的流动是一类复杂流动,并在石油开发、能源安全以及地下水污染治理等许多领域均有涉及。在这些领域中,很多问题需要根据已有经验和已测数据构建地下流体模型,并对地下流体的某些属性或发展趋势做出预测,例如水油比、CO浓度、放射性粒子浓度或者生产曲线等。在石油储层的勘探开发过程中,地下流量的预测尤为重要,不仅关乎生产效益更对生产决策起到至关重要的作用。

在地质流体力学中,以孔隙度和渗透率作为参数的流量概率密度分布,往往呈现出较为复杂的函数结构,很难用常用的经典分布(如标准正态分布,均匀分布等)予以描述。因此,需要从孔隙度、渗透率及油井产量所构成的复杂参数/分布结构中获得必要的模拟样本。该文将马尔可夫链蒙特卡罗方法(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方法的基本概念及相关算法,并依据国外相关研究实例介绍了其在多孔介质中流体预测中的运用。详细阐述了通过MCMC方法获得多孔介质流体模型中渗透率和孔隙度的贝叶斯估计的基本步骤,进而依据其估计建立出完整的地下属性描述并进行流量预测的具体过程。同时在文章的结尾部分,根据现有方法的不足,提出了相关的改进设想,有望基于相关分析提高预测成果。

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

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方法篇4

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方法篇5

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

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

MCMC方法篇6

股市的金融计量模型是从研究股市波动率开始的, Engle (1982) 提出ARCH模型[1], 对波动率聚类的进行了建模, 这个结果迅速后被用于金融时间序列的分析, 并提出一系列的模型, 形成GARCH模型族。

GARCH模型族面临的一个质疑是它没有考虑金融时间序列因为结构突变所形成的不同状态[2]。Lamoureux等 (1990) 对收益率数据加入结构突变, 进行模拟试验后得出, 结构突变能显著降低GARCH模型的持久性参数[3]。因而, 不能排除波动的聚类还来源于结构突变。Smith (2002) 和Maekawa等 (2005) 分别在美国和日本股市中发现明显的结构突变现象[4,5]。

之后, Hamilton和Susmel (1994) [6], Cai (1994) [7]最先结合结构突变, 将马尔可夫转换模型运用于股市波动性分析, 提出马尔可夫转换-ARCH模型 (Markov Switching ARCH Model) 。观察到股市波动率高的状态往往是经济比较低迷时, Hamilton和Lin (1996) 又结合经济周期, 用双变量马尔可夫转换模型进一步估计股市波动性, 得出经济萧条是股市高波动率的最重要的影响因素[8]。在亚洲股市的研究中, Fong (1997) 和Li和Lin (2003) 分别分析了日本股市和台湾股市的波动性[9,10]。

上述马尔可夫转换模型都是用极大似然法估计的, 但是极大似然法有很大的局限, 比如它不能确定是否找到了最高似然值。而另一种估计方法——马尔可夫链蒙特卡罗法 (MCMC Method) , 可以很好地弥补这个缺陷。自从Albert和Chib (1993) 将MCMC模拟法首度用于马尔可夫转换模型估计[11]以来, 该方法的应用发展很快。Kaufmann和Fruehwirth-Schnatter (2002) 用此方法估计了马尔可夫转换ARCH模型[12];Smith (2002) 估计了马尔可夫转换随机波动模型[13]。Das和Yoo (2004) , Henneke等 (2006) 分别估计了MS-GARCH模型[14,15]。

研究中国股市的金融计量模型中, 孙金丽等 (2003) 、蒋祥林等 (2004) 、万军等 (2008) 、江孝感等 (2009) 都用极大似然法估计了马尔可夫转换-GARCH模型, 他们都说明了此前学术界所发现的GARCH模型持久性参数过大的问题, 其中前两者发现中国股市的波动性状态转换同政策相关, 从而提出中国股市的波动主要由政策引起[16,17,18,19]。本文将运用MCMC方法估计上证综指的马尔可夫转换-ARCH模型, 并且得出更为丰富的结论。

2 马尔可夫转换-ARCH模型

本文模型参考Henneke等 (2006) 和Rachev等 (2008) , 假设n个状态, 每个时间点所处状态St随机出现, 状态之间的转换服从马尔可夫转移概率矩阵, 记πij是从状态i到状态j的转移概率。马尔可夫转换-ARCH模型表述如下:

$\begin{array}{l} y_{t} = ω_{S_{t}} + ε_{t} (1) \\ h_{t} = α_{S_{t}} + \sum_{q = 1}^{Q} β_{q, S_{t - q}} ε_{t - q}^{2} (2) \end{array}$

其中, εt～N (0, ht) , ωSt, αSt, {βq, St}Qq=1各参数依赖于t时期的状态St, ht是条件方差序列。

一般而言, 方程 (1) 又称为均值方程, 而方程 (2) 称为方差方程。其中均值方程没有包括任何自回归项或移动平均项, 因为在对AR-ARCH建模过程中, 本文发现εt服从t分布的假设将使自回归项不再具有统计显著性。

3 模型的MCMC估计

3.1 MCMC估计方法

首先设定各组参数的先验概率密度函数。对于ARCH参数, 通常采用正态分布假设,

$θ_{S} \sim Ν (μ_{S}, V_{S}) Ι_{{θ_{S}}}$

其中, I{θS}是GARCH参数的示性函数, 在满足GARCH参数的限制时, 该示性函数等于1, 否则等于0, 即

$Ι_{{θ_{S}}} = {\begin{cases} 1, α_{S} > 0, {β_{q, S}}_{q = 1}^{Q} \geq 0, α_{S} + \sum_{q = 1}^{Q} β_{q, S} < 1 \\ 0, α_{S} \leq 0, {β_{q, S}}_{q = 1}^{Q} < 0, α_{S} + \sum_{q = 1}^{Q} β_{q, S} < 1 \end{cases} (3)$

转移概率πi的先验概率密度函数设定为

$π_{i} \sim D i r i c h l e t (a_{i 1}, a_{i 2}, \dots, a_{i S}) (4)$

为了进行MCMC模拟, 先求出后验概率密度函数, 然后分别计算状态序列St和转移概率πi的条件概率密度函数, 以及GARCH参数的似然函数。

后验概率密度函数

$p (θ, π, S | y) = p (y | θ, π, S) p (θ) p (π) (5)$

其自然对数函数为

$\begin{array}{l} l o g p (θ, π, S | y) \\ = c - \frac{1}{2} \sum_{t = q + 1}^{Τ} [l o g h_{t} + \frac{(y_{t} - ω_{S_{t}})^{2}}{h_{t}}] \\ - \frac{1}{2} (θ_{S} - μ_{S})^{'} V_{S}^{- 1} (θ_{S} - μ_{S}) Ι_{(θ_{S})} \\ + \sum_{i = 1}^{S} \sum_{j = 1}^{S} (a_{i j} + n_{i j} - 1) \log (π_{i j}) \end{array} (6)$

其中, c为常数项, nij为状态序列中相邻时间点从状态i向状态j转换的次数。式 (6) 显示转移概率向量πi只同状态序列{St}Tt=1有关, 即

$\log (π_{i} | y, Φ_{- π_{i}}) = c_{π} + \sum_{j = 1}^{S} (a_{i j} + n_{i j} - 1) \log (π_{i j}) (7)$

其中, cπ是不重要的常数, Φ-πi指模型除了πi的所有参数。这个条件概率密度函数就是Dirichlet分布函数的核的对数值。

根据贝叶斯法则, St的条件概率密度函数

$\begin{array}{l} p (S_{t} | y, Φ_{- S}, S_{- t}) \\ = \frac{p (S_{t} = i, S_{- t}, y | Φ_{- S})}{p (S_{- t}, y | Φ_{- S})} \\ = \frac{p (y | Φ_{- S}, S_{- t}, S_{t} = i) p (S_{t} = i, S_{- t} | Φ_{- S})}{p (S_{- t}, y | Φ_{- S})} \\ = \frac{p (y | Φ_{- S}, S_{- t}, S_{t} = i) π_{j i} π_{i k}}{\sum_{s = 1}^{S} p (y | Φ_{- S}, S_{- t}, S_{t} = s) π_{j s} π_{s k}} \end{array} (8)$

ARCH参数的后验概率分布并不是标准形式, 其核的对数函数为

$\begin{array}{l} \log (p (θ_{S_{t} = i} | Φ_{- (S_{t} = i)}, y)) \\ = c_{θ} - \frac{1}{2} \sum_{t = q + 1}^{Τ} [l o g h_{S_{t} = i} + \frac{(y_{t} - ω_{S_{t} = i})^{2}}{h_{S_{t} = i}}] \\ - \frac{1}{2} \sum_{s = 1}^{S} (θ_{s} - μ_{s})^{'} V_{S}^{- 1} (θ_{S} - μ_{S}) Ι_{(S_{t} = s)} \end{array} (9)$

因为ωSt通过残差影响条件方差ht, ωSt将出现在方括号中的三个地方, 使式 (9) 复杂化。ARCH参数的条件概率分布不再是标准分布, 为了得到该后验概率分布的随机取样值, 本文用随机游走Metropolis-Hasting取样法。

MCMC模拟将MS-ARCH模型的估计分解成以下几步:

第一步:基于ARCH参数θ (j) St和转移概率π (j) i, 通过状态序列St的条件概率密度函数取得第 (j+1) 次随机序列{St} $_{t = 1}^{Τ}$ (j+1) 的模拟值;

第二步:基于θ (j) St和{St} $_{t = 1}^{Τ}$ (j+1) , 通过转移概率πi的条件概率密度函数取得转移概率的第 (j+1) 次模拟值π (j+1) i;

第三步:基于{St} $_{t = 1}^{Τ}$ (j+1) 和π (j+1) i, 通过Metropolis-Hasting取样法得到ARCH参数的第 (j+1) 次模拟值θ (j+1) St.

不断重复以上三步各参数的模拟, 直至参数模拟值各自收敛, 这些收敛之后的参数模拟数的平均值和标准差就是MCMC模拟法的参数估计结果。

3.2 模型的MCMC估计结果

本文上证综指的数据来源于Wind资讯, 共采集了1994年12月30日至2009年10月15日的共740个周收盘价, 然后通过收盘价的自然对数之差得到739个连续周收益率。为了消除极值的负面影响, 本文把这些数据中绝对值大于0.15的收益率都变为0.15, 这样的极值共有4个, 其中三个发生在1995年, 一个发生在1996年。1996年12月之后, 由于中国证监局实行了涨跌停板制, 就再没有发生周涨跌幅超过15%的波动了。

根据经验, MS-ARCH模型的状态数一般较小, 不会超过3个。因此, 本文只估计了两状态和三状态MS-ARCH模型。本文用Matlab软件编程, 共进行了30000次模拟。对模拟结果进行Smith和Robert (1992) 的收敛性检验证实参数收敛。至于参数随机值的接受概率, 第一、第二和第三状态分别为22%、48%和46%。三个接受概率都在20%到75%的最佳概率区间。

将参数随机中前面10000个模拟值作为废弃数据 (burn-in data) 处理, 并根据后面20000个模拟数据求参数的平均值和标准差, 得到的结果如表1所示。为了比较模型, 本文同时估计了ARCH (3) 模型。

表1显示了三个模型的估计结果。各个状态的α项具有5%的统计显著性, 其他参数除了αS=2皆不显著。但是, 每个状态的β参数相加则具有统计显著性, 比如MS (3) -ARCH (3) 中βS=1均值为0.3916, 标准差为0.1617。这样, 虽然β参数没有很好的稳定性, 但每个状态的无条件方差却非常稳定。

注: 括号内数值是参数的标准差。**、*分别表示在5%和10%水平上统计显著性, 下同。

3.3 MS-ARCH模型的比较

表1显示, 三个状态的MS-ARCH模型比其他两个模型具有更高的似然值、AIC和BIC信息准则, 因此更优越。另外, 本文估计了GARCH族模型中主要模型, 并求出各自的似然值和AIC, BIC信息值, 估计结果如表2。表2中的四个GARCH族模型比较简洁, 并且似然值在同类中比较高, 各参数也基本具有统计显著性, 因此这四个模型代表了GARCH族中最好的模型。

然而, 在同MS-ARCH模型的比较中, 似然比检验也显示, 考虑了结构突变的后者明显优于这些普通的GARCH族模型。故而, 按照AIC和BIC信息准则, 三状态MS-ARCH模型优于几乎所有的普通GARCH族模型。

3.4 上证综指MS-ARCH模型的分析

表1显示, 在考虑了结构突变之后GARCH模型的持久性指标降低, 显示该指标难以有效说明波动持续, 而状态的持久性才是造成股市波动聚类的主要原因。波动状态的持续性对波动聚类有更强的解释能力, 并且不同状态的波动聚类性质也不一样。

三个状态的无条件方差分别为0.0013、0.0025、0.0049。其中, 第三个状态的无条件方差是第一个状态的近4倍。为此, 第一、第二和第三状态又分别称为低波动状态、中等波动状态和高波动状态。另外, ARCH (3) 模型的期望条件方差为0.0028, 介于低波动和高波动状态之间。通过模型所估计出的转移概率矩阵, 本文发现, 低波动状态的持续概率最大, 而高波动状态最不易持续。其中, 低、中等和高波动三个状态的平均持续时间分别为21.3周、8.3周和6.3周。

通过比较三个状态在每个时间点的各自发生概率, 本文观察到, 低波动状态平均持续时期最长, 并且发生的总时间也最长, 而高波动状态持续时期短, 也不经常发生。比较有意思的是, 中等波动状态基本伴随着牛市而来, 如1997年初、1999年上半年、2005年股改之后的大牛市以及之后2009年的反弹。这一点从表1的估计结果中也可以得到印证, 即ωS=2显著大于0。另外, 本文发现高波动性同熊市以及股市拐点处的市场相关性非常大。

下面比较ARCH (3) 和三状态MS-ARCH模型的条件方差ht, ARCH. 如图1和图2显示, 三状态MS-ARCH模型的条件方差序列ht, MS中低波动状态的方差较小, 而中等和高波动状态的方差普遍大于ht, ARCH的平均方差, 这一点表1的估计结果也可证实。然而, ht, MS最明显之处在于它的方差往往有聚类倾向, 这是中等和高波动状态发生时, 波动率保持高位的表现。对照图2和中等和高波动状态发生的概率, 方差聚类发生的时期, 同中等或高波动状态的概率接近于1的时期高度相关。

这两张条件方差图显示, MS-ARCH模型把方差按不同的水平分成几个层级, 每个层级也相当于一个波动状态, 然后在每个状态分别进行ARCH模型估计, 而让状态之间的转换服从马尔可夫链式的转移概率。总体而言, 这种方法能显著提高波动率的拟合程度。

本文将上证综指的MS-ARCH估计结果同世界其他股市的相似模型进行比较, 并且发现一些有趣的结论。本文模型中, 均值方程的自回归项系数没有统计显著性, 因而本文在建模时没有加进自回归项。然而, 台湾、日本和美国股市却表现出一阶自回归系数的统计显著性。至于方差方程中β参数的没有统计显著性, 在其他股市如台湾、日本和美国也有相似的情况。

模型结果显示, 上海股市的高波动状态波动率是低波动状态的将近4倍。而该数字在台湾股市是15.62倍, 美国是13倍, 在日本股市中两状态模型的波动率相差6倍。然而, 上海股市的4倍并不是说明中国大陆股市的高波动状态的波动率不如其他股市大, 而是, 上海股市的低波动状态的波动绝对值也相对较大, 显得差距倍数较小, 而这正是中国大陆股市的总体波动性显著较大的一个证据。

在高波动状态的持续时间的比较上, 上海股市相比台湾、日本以及美国的股市, 我国股市的高波动状态持续时间却显得短。台湾地区的低、中等和高波动状态持续时间分别为114周、72周和60周, 明显比上海股市长。美国分别为132周、116周和99周。日本低波动和高波动状态的平均持续时间分别为40周和22周。由此, 本文得出, 上海股市相对于其他股市波动状态转换频繁, 高波动状态的持续周期竟然只有6周。究其原因, 上海股市的波动是投机者心理波动的反应, 而不是周边政治经济环境波动的反应。

同中国大陆股市的中等和高波动状态分别与牛熊市相关不一样, 台湾的高波动状态却同台湾地区以及周边的政治和经济危机有关, 如1996年大陆在台海的导弹发射、1997年的亚洲金融危机、1998年资本所得税的征收等。日本的高波动状态也同国际政治经济危机相关, 他的两个主要高波动状态分别是1986～1988年和1990～1993年。前一时期因为广场协议之后日元急剧升值引起日本经济下滑, 再加上1987年的华尔街股灾使日本的股市波动率加大。后一时期中, 海湾石油危机和日本股市泡沫破裂后所造成的股市持续的高波动性。美国的高波动状态同经济萧条期息息相关, 比如1962年下半年、1969年5月至1972年1月、1973年2月至1976年4月, 以及1987年股灾后。在这些国家和地区的股市都发挥着政治和经济晴雨表的功能, 而中国大陆股市却没有。

4 结论和建议

本文通过MCMC模拟法估计了上证综指MS-ARCH模型。估计结果显示, 三状态MS-ARCH模型好于两状态模型, 并且按照模型似然值和AIC, BIC信息准则, 本文的模型优于几乎所有的普通GARCH族模型。

除此之外, 本文模型更深刻的解释了上海股市波动的持续性, 波幅, 转换频率以及形成原因, 然后又将上海股市的相关结论同世界其他地方的股市相比较, 得出富有启发性的见解。本文的主要结论如下:

①波动率的持续性。上海股市各个波动状态的平均持续时间都偏短, 相比其他股市的低波动状态能持续数年, 上海股市表现出频繁的波动状态转换。

②波动率的大小。模型显示上海股市整体而言, 比台湾、日本和美国的波动率大, 但差距倍数不如其他股市大。

③波动状态形成的原因。台湾股市的高波动状态反映了台湾地区的政治经济危机, 日本股市的高波动状态同日本的世界经济状况息息相关, 而美国的高波动状态往往伴随着美国经济的萧条期。而本文发现上海股市的中等波动状态同上证综指的牛市相关性较明显, 而高波动性往往伴随着熊市或调整市。这都说明, 上海股市并不是中国经济的晴雨表。

④收益率和波动率的关系。中等波动状态的平均周收益率达到2.11%, 具有明显的统计显著性。中国股市中, 收益率和波动率在中等波动状态居然有如此强的相关关系, 是本文最大的发现。这个结果显示, 中等波动到来的时间点具有很强的实际应用价值, 可用于构建数量投资策略。

至于政策建议, 本文提倡科学的数量投资法, 使市场逐步减少投机成分, 增加理性科学的投资方法, 从而逐步促进股市的完善。

摘要：引入结构突变, 对上证综指马尔可夫转换-ARCH模型通过马尔可夫蒙特卡罗方法 (MCMC方法) 进行估计。在30000次参数模拟之后, 本文得到稳健、可靠的结果, 似然比检验显示本文模型好于几乎所有GARCH族模型。本文结论: (1) 相对于世界主要股市, 中国股市各波动状态的持续时间短、波动幅度大; (2) 不像其他股市, 中国股市的波动不能反应国内外的政治经济状况的变化; (3) 中国股市中等波动状态的收益率显著大于0。这些结论提供了一个认识中国股市波动性的全新视角, 还揭示了一种基于模型的实用数量投资方法, 最后本文提出了完善中国股市的相关建议。

MCMC方法篇7

复杂的接收机检测算法是影响MIMO系统大规模商用的一个原因。理论上,最大后验概率MAP(Maximum A Posterior)检测算法是最优的检测方法,但因其复杂度太高而很难实际应用[2]。MIMO系统中利用软信息的迭代检测算法因其优良的性能而成为研究的热点。

近年来,MCMC算法因其较低的复杂度而被应用在MIMO系统中[3,4,5,6]。参考文献[6]设计了一种基于MCMC算法的迭代接收机,利用贝叶斯判决这一统计方法来获得低误码率和高性能。本文提出一种融合了SIC-MMSE算法的MCMC算法,可以有效利用反馈软信息,仿真结果表明,该算法可以有效降低传统MCMC算法的复杂度。

1系统框图

图1描述了MIMO系统的带迭代检测和译码的接收机结构。假设一个发送天线数为M接收天线数为N的MIMO平坦快衰落信道,系统模型可表示为:

其中,s=[s1,s2,…sM]是发送信号符号向量,其元素属于星座点集合S=[α1,…αj,…αMc],其中αj(1燮j燮Mc)为数字调制的星座点,Mc为星座点集合S的大小,表示每个符号对应映射的比特数,4-QAM时,Mc=2,16-QAM时,Mc=4。r是N×1维接收信号向量,H是N×M维信道传输矩阵,其元素Hi,j(i∈[1,N],j∈[1,M])表示从发送天线j到接收天线i的复路径增益。n是N维零均值高斯噪声向量,其协方差矩阵为E[nn H]=σn2INr,IN为N×N维单位阵。

每个码元间隔内发送的信息向量b中包含McNt个信息比特,记为b=[b1,b2,…bMc Nt]T。在讨论MIMO软判决检测时,通常使用+1和-1来分别表示逻辑1和逻辑0。对于这样的比特向量用x=[x1,x2,…xMcNt]T来表示。最大化一个比特位的最大先验概率或者先验概率可以使此位判断错误的概率最小,先验概率通常表达为对数似然值(LLR)ln P(xk=+1)/P(xk=-1),通过观察LLR的符号可以判决该比特值,LLR的绝对值体现了此次判决的可信度。简单的加减运算即允许从在译码循环中得到的新(外)信息里分离出先验信息或称旧信息。

传统的接收机结构由SISO(软输入软输出)检测器和FEC(前向纠错)信道译码器两个基本处理器组成,两个检测器之间交换软信息[2]。本文利用SIC-MMSE算法的辅助,改变了吉布斯采样传统的随机采样方式。基于接收信号序列r,上次迭代由信道译码器传回的先验信息LAb以及预判决器输出的发送信号预估计值x赞,检测器生成符号序列的软输出序列。在检测器的输出LDb中减去先验信息后,相对信道译码器是新(外)信息的剩余信息LEb解交织,成为信道译码器的先验输入LAc。类似地,在信道译码器的输出里分离出其软输入信息,从而生成新(外)信息,经交织传送给预判决器,作为先验信息LAb和发送信号预估计值x赞一起反馈给检测器。

2算法描述

定义向量x-k=[x1,…xk-1,xk+1,…,xMc Nt],LEc=[LEc(x1),LEc(x2),…LEc(XMc Nt)]即每个码元间隔的全部发送比特的外信息。已知接收向量r,则:

其中,,需要计算所有可能的x-k。而x-k可能的组合个数按McNt的指数增长。这个复杂度问题可用蒙特卡洛综合方法[3]来解决。

2.1蒙特卡洛综合

考虑一个估计随机变量X的函数h(x)的加权均值的一般性问题,给定加权函数f(x),亦即:

其中χ是X的定义域,f(·)是密度函数,即f(x)叟0,并且。由参考文献[3],可以通过对经验均值的估计来估计上式,得到:

这种方法的收敛速度可以由h的方差来估计:

可以推断出其估计精度随着采样点数的平方下降而下降。因此,为了提高估计精度通常采用增加采样点个数的方式。在其他条件不变的情况下,采样点个数的增加的确可以带来估计精度的提高。而本文提出了一种不增加采样点数,而是改变采样点初值的方法使得式(5的分子变小,进而在不增加采样点数甚至有可能减少采样点数的情况下,提高式(4)的近似程度。

2.2吉布斯采样

本文中应用MCMC方法由P(x|r,LEc)来产生采样点。K(x,x′)定义了由x指定的状态到由x′指定的状态的转移。马尔科夫链中的状态数按x的大小指数增长,吉布斯采样[4]是MCMC检测器产生发送信息序列x的采样的一般方法,用来缓解复杂度问题。预判决器的输出和先验信息送入SISO检测器,经Ns次迭代运算产生Ns组采样序列,用来计算LLR(对数似然概率)值。本文中,检测器采用的一种吉布斯采样步骤如下:

将(7)代入(6),计算出。以此概率抽取xk(n)的方法[5]是:随机抽取一个0和1之间均匀分布的数u,若u小于,则xk(n)取+1,否则xk(n)取-1。

2.3 LLR的计算

由参考文献[3]知外信息较精确的计算公式为:

经简化得LEb(xk)的近似计算式[6,7]为:

计算出的LEb(xk)传递给信道译码器。从(4)式也可以直观地看出,采样点数目对整个算法的复杂度有很大的影响,因此若在吉布斯采样时减少采样点可以使得后来的运算量减少。

2.4 SIC-MMSE原理

本文对传统方法的改进是在吉布斯采样的过程中,借助了SIC-MMSE算法[8,9],从而使得采样点更加有效,马尔科夫链迅速收敛。得到由信道译码器返回的软信息后,将其送入预判决器,根据得到的LLR值,可以定义第m个发送天线的相关软输出符号,更精确的软输出符号的均值为:

其中,q是调制量座点数,对M-QAM,q=M,sm(q)表示符号sm的第q个合理值。

对于16-QAM调制方式,有:

对于常系数调制方案,定义第m个天线发送的符号的协方差为:

根据软干扰消除原则,由MMSE算法计算得到的对第m个发送天线的发送符号的估计可以表达为:

其中,MMSE加权矩阵WMMSE的第m列可以表示为:

其中,Ip代表P行P列的单位矩阵,V=diag[v1,v2,…,vM]。

融入SIC-MMSE算法后,由LLR经过预判决得到s赞m的初始化x(0),从而使得马尔科夫链收敛加快,减少采样点,提高估计精度,降低算法复杂度。

3仿真结果及分析

仿真采用天线数为4×4(4个发送天线,4个接收天线)的MIMO系统,信道编码为码率R=1/3的Turbo编码。仿真设定MIMO检测器为4次迭代,Turbo译码器内部为8次迭代,信息序列长度为9 216 bit。图2和图3分别比较了16 QAM和64 QAM调制方式下传统的MCMC算法和改进的MCMC算法的误码率。

图2和图3证明了2.1节的推断,即采样点数的增加可以带来MCMC算法精度的提高,从而降低误码率,并且证明了本文算法取得了相对传统算法更优的性能。如图2中,取Ns=5,当误码率达到10-5时,本文算法比传统的MCMC算法性能提高了约0.5 dB,图3中,取Ns=10,当误码率达到10-5时,本文算法性能提高了约2.5 dB。另一方面,从图2和图3的仿真结果可以得出要取得相同的性能,本文算法所需的采样点少于传统的MCMC算法的结论。例如,要实现4 dB之前达到目标误码率,传统的MCMC算法需取Ns=20,而本文算法取Ns=10即可。

表1列出了SIC-MMSE算法辅助的MCMC算法和传统算法在帧长1 024、16 QAM调制的情况下,取不同采样点数时,达到误码率为10-5时的SNR点和总运算次数。

由表1可以看出,由SIC-MMSE辅助的算法可以以相同的复杂度取得更高的性能。例如,由SIC-MMSE辅助的MCMC算法取Ns=10,而传统的MCMC算法取Ns=20时,都实现了4 dB之前达到目标误码率,而此时前者的复杂度近似为后者的二分之一。

本文提出了一种MIMO系统中低复杂度的MCMC检测算法,该算法在传统的MCMC算法基础上,融合了SIC-MMSE算法,更加有效利用信道译码器产生的软信息,提高蒙特卡洛综合方法的近似程度,可以在不损失检测性能的前提下减少采样点数,从而降低了MIMO接收机的复杂度。仿真结果显示,由SIC-MMSE辅助的MCMC算法与传统的MCMC算法比较,在取得相同的误码率时,系统复杂度低很多。综合考虑复杂度和检测性能,由SIC-MMSE辅助的MCMC算法明显优于传统的MCMC算法。

参考文献

[1]FOSCHINI G J.Layered space-time architecture for wire-less communication in a fading environment when using multi-element antennas[J].Bell Labs.Tech.j.,1996,1(2):41-59.

[2]HOCHWALD B M,BRINK S T.Achieving near-capacity on a multiple-antenna channel[J].IEEE Transactions on Communications,2003,51:389-399.

[3]BEHROUZ F B,ZHU Hai Dong,SHI Zhen Ning.Markov chain monte carlo algorithms for CDMA and MIMO communication systems[J].IEEE Transactions of signalprocessing,2006,54(5):1896-1906.

[4]HENRIKSEN S,NINNESS B,WELLER S R.Convergence of Markov-Chain Monte-Carlo approaches to multiuser and MIMO detection[J].Selected Areas in Communica-tions,IEEE Journal on,2008,26(3):497-505.

[5]ZHU Hai Dong,BOROUJENY B F.MIMO Detection using markov chain monte carlo techniques for near-capacity performance[J].Acoustics,Speech,and Signal Processing,2005.Proceedings.(ICASSP'05).IEEE Inter-national Conference on,2005,3:1017-1020.

[6]ZHU Hai Dong,BOROUJENY B F,CHEN Rong Rong.On performance of sphere decoding and markov chain monte carlo detection methods[J].Signal Processing Letters,IEEE,2005,12(10):669-672.

[7]LARAWAY S A,BOROUJENY B F.Implementation of a markov chain monte carlo based multiuser/MIMO detector[J].Circuits and Systems I:Regular Papers,IEEE Trans-actions on,2009,56(1):246-255.

[8]TUCHLER M,SINGER A C,KOETTER R.Minimum mean squared error equalization using a priori information[J].Signal Processing,IEEE Transactions on Acoustics,Speech,and Signal Processing,2002,50:673-683.

MCMC方法篇8

1 计算机图像处理技术和图像处理过程的解析

计算机图像处理技术作为推动现代化信息技术发展的核心技术手段, 在实际应用过程中是非常常见的。简单来说, 计算机图像处理技术其实就是一个二维数据矩阵生成的一个过程。在一般情况下, 计算机当中的二维矩阵通常想要实现构成, 那么必须要经过采集、扫描和量化这三个环节。其中, 当计算机在进行扫描作业时, 必须要注重计算机对图像的扫描序列制定, 必须要和事先制定的相关图像序列保持绝对的一致。另外, 计算机在处理图像二维矩阵过程中, 主要就是对图像像素当中的灰度值进行准确的测量, 而像素值则是扫描过程中最小的数据单位, 在一般的情况下, 计算机的图像处理步骤通常都是借助光电传感技术完成的。这里我们可以注意, 计算机图像获取的过程其实就是对图像灰度值的获取, 而且最终将图像的灰度值转换为图像的理算数值。所以, 计算机图像处理的过程就可以看做是在一个矩阵空间内对图像网格进行扫描处理, 进而获得一个相应的整数矩阵, 而这个整数举证与图像扫描之间的关系是相互交错的, 因此也可以说是计算机图像是对于图像数字化之后产生的二维图像矩阵数据。因此, 计算机图像处理的整个过程就为, 景物确定→成像处理→图像分析→图像采样→图像量化处理→数字图像。

2 Bayesian-MCMC算法在计算机图像处理中的应用措施

2.1 Bayesian算法在计算机图像处理流程中的应用

Bayesian算法又可以被称为是贝叶斯分类计算处理技术, 主要应用于计算机图像处理过程当中, 技术的运转主要是依据结合叶背斯定理实现的, 其计算公式如 (1) 所示, 这就是叶背斯分类计算公式的主要特点。

当计算机在英语Bayesian-MCMC算法对图像进行应用处理时, 其主要的方式就是按照特定的图像计算规律对其进行随机分布数据矩阵, 其中每一个单元参数都是随机变化量, 并且这个变化量不属于某一个特定的变化类别, 同时也不属于多个不相容的参量变化事件, 所以, 当计算机在处理蚃信息时, 就可以使用Bayesian算法来对其进行相应的计算处理, 因此, 在这种计算处理情况之下, 我们就可以将公式 (1) 转换为公式 (2) , 如下所示:

从变化得出的公式 (2) 看到, 其中数值i可以分别表示为1、2、3、4.....N, 而公式当中的X则是图像信息在抽取相应特点数值和信息分类之后得出的图像灰度值, 数据P则充分表明了图像当中的概率条件信息, 因此, 在得知了图像归属概率参数的情况之下, 就可以根据这一条件系数进行分类计算处理, 最终达到对图像信息进行计算的目的, 并且满足相应的图像计算要求。

2.2 MCMC算法在计算机图像处理中的应用

MCMC算法是结合了随机模拟处理技术和马尔卡夫链原理的一种计算应用技术。依据贝叶斯分类计算的相关技术原则, 在对图像进行实际密度计算之后, 可以根据图像矩阵数据参数、正则因子等, 进而对其进行随机模拟计算处理。

2.3 Bayesian-MCMC算法在计算机图像中的应用

根据Bayesian算法和MCMC算法在计算机图像中的应用特点来看, 综合使用两种图像处理技术能够大大改善单一图像处理过程中的缺陷, 其图像处理步骤如下所示:

(1) 首先, 必须在某一时刻找到合适的马尔卡夫链条并创建连接, 同时确定数据核的转移, 之后确定某一区域的平稳分布状况。 (2) 随后从马尔卡夫链中的某一节点出发, 然后再通过上述步骤来建立马尔卡夫链的计算处理序列。 (3) 最后在实际计算处理过程中, 应该注意图像预算迭代数值的处理, 根据这一图像处理过程就可以得出图像预算公式 (3) , 公式 (3) 也可以看做是Bayesian-MCMC算法应用计算机图像处理计算的原则。

3 结语

总之, 在计算机图像处理过程中应用Bayesian-MCMC算法对图像计算处理的促进性是显而易见的, 不仅能够大大降低图像处理计算的整个过程的难度, 还可以让图像计算处理具有较高的效率, 其实际应用的优势非常明显。

参考文献

[1]杨海东, 肖宜, 王卓民, 邵东国, 刘碧玉.突发性水污染事件溯源方法[J].水科学进展, 2014 (01) :122-129.

[2]项立.Bayesian-MCMC算法在计算机图像处理中的应用[J].河南科技, 2014 (06) :13-14.

【MCMC方法】推荐阅读：