马尔可夫（精选8篇）
马尔可夫 篇1
脑部磁共振 (magnetic resonance, MR) 成像技术广泛应用于临床实践和神经科学研究领域, 其能提供良好的空间分辨率和软组织对比度[1]。将脑部磁共振 (MR) 图像分割出三个重要组织:灰质 (GM) 、白质 (WM) 和脑脊液 (CSF) , 这些组织的组合物在指定区域内的变化可用于表征生理过程和疾病严重程度, 由此可见, 脑部MR图像分割具有重要的研究意义。
高斯混合模型 (GMM) [2]是神经学中广泛使用的一种统计学模型, GMM不仅与脑部MR图像的分段常数性质一致, 而且具有较低的计算复杂度。GMM的模型参数可以利用期望最大化 (EM) 算法根据最大似然 (ML) 准则来估计, 然而, 基于EM的ML估计具有过度拟合和容易限于局部最优解的缺点。为了克服这些缺点, 使用几种全局优化技术来替代EM算法, 例如, 文献[3]在似然估计中结合了遗传算法, 提出了GA-EM算法。此外, 当先验知识可用时, 最大后验概率 (MAP) 估计是ML估计的一种常见替代方法。文献[4]提出一种MAP-MRF框架来求解图像分割问题, 通过将体素类标签建模为马尔可夫随机场 (MRF) 来表示体素空间依赖性的先验。文献[5]提出一种基于全局随机搜索的推理方法, 即马尔可夫链蒙特卡尔 (Markov chain monte carl, MCMC) 推理, 用来替代确定性程序。文献[6]受免疫机制启发, 提出一种克隆选择算法 (clonal selection algorithm, CSA) , 基于克隆选择理论, 选择能够识别抗原的抗体来进行繁殖, 繁殖的细胞会通过一个亲和力成熟过程来改进它们对抗原的亲和力, CSA模仿对抗原刺激免疫应答机制来实现全局最优。
将CSA和MCMC技术融合到隐马尔可夫随机场 (hidden Markov random field, HMRF) 模型估计中, 提出一种用于脑部MR图像分割的HMRF-CSA算法。首先, 通过MCMC方法近似最优标签配置, 然后, 由CSA算法估计HMRF模型参数。用全局随机优化技术替代确定性搜索程序, 以此提高分割算法的鲁棒性, 同时, 将MR图像建模为分段常数图像的乘法分量, 根据MCMC推断方法获得的中间分割结果来评估图像的不均匀性。通过仿真实验, 将本文HMRF-CSA算法与现有的GA-EM方法、可变形共同分割 (D-C) 算法、SPM软件包中的统一分割算法和FMRIB软件库 (FSL) 上的HMRF-EM分割算法进行比较, 结果表明该算法具有更好的分割精度。
1 相关技术
1.1 图像不均匀模式
由图像采集的不完善所导致的图像不均匀性, 或称偏场或强度非均匀性 (INU) [7]是MR图像分析的难点之一。设定y={yi;i=1, 2, …, N}表示一副脑部MR图像, 其中yi表示在体素i处的强度, N表示体素的数目, 未知偏场B={bi;i=1, 2, …, N}通常建模为y的乘法分量, 如下式所示
式 (1) 中, 是理想图像, 是附加的高斯白噪声。由于图像中偏场B变化缓慢, 所以可将它定义为在整个图像域上的一个平滑函数。采用正交多项式{Wj:j=1, 2, …, NOP}作为偏置函数来近似偏场[8]
式 (2) 中, φ={φj:j=1, 2, …, NOP}表示组合系数, NOP= (D+1) (D+2) /2是多项式的数目, D是多项式的度。
1.2 统计学模型
假设体素强度y={yj;j=1, 2, …, N}符合GMM;从有先验概率πk的高斯分布N (μkΣk) 中独立采样每个强度值yj, 观察图像的似然, 计算式如下:
通过最大化上述似然函数来估计最优GMM参数, 确定参数后, 利用贝叶斯分类器对每个体素进行分类, 以此求解脑部图像分割问题[7]。
为了将空间约束融入到这个模型中, 本文应用MRF到模型类标签x={xj;j=1, 2, …, N}中, 根据Hammersley-Clifford理论, 类标签p (x) 的先验联合分布符合Gibbs分布。在MAP-MRF框架下, 图像分割等价于通过最大化其后验概率寻找最优配置x*
式 (4) 中, Θ={μk, Σk;k=1, 2, …, K}表示模型参数, p (y|x;Θ) 是图像似然, p (x) 是空间先验。
本文将图像强度y当作在相同图像点阵中另一个随机场建模的模型, 然后将代表潜类标签的MRFx变成HMRF。在这个模型中, 将图像分割问题制定为配置x和参数Θ的最大联合概率
式 (4) 中后验与式 (5) 中联合概率之间的差是惩罚项p (Θ|y, x) , 用于检查模型参数是否与配置x给出的观察值一致。
2 HMRF-CSA算法提出
HMRF模型主要用来估计式 (5) 中的最优类标签和模型参数 (x, Θ) , 估计过程可划分成两个相互依赖的优化步骤:搜索最优配置x*和学习最匹配模型参数Θ*。使用下式三个步骤的迭代程序来实现HMRF模型估计
式 (6) 中, f (., .) 是基于观察y和分割结果x纠正偏场的函数, t∈{1, 2, …, Tmax}表示当前迭代数目。在每次迭代中, 首先采用MCMC方法实现MRF-MAP估计, 在近似的分割结果下估计偏场, 然后利用CSA学习HMRF模型参数, 当达到最大迭代数目或分割结果变成稳态时迭代停止。
2.1 MCMC体素分类
上述迭代步骤中第一步是通过MRF-MAP近似寻找最优配置x*, 使用MCMC方法求解这种优化问题, 根据式 (4) , 对于给定任意特定配置x, 假设yj相互独立且符合基于参数Θk={μk, Σk}的多元高斯分布, 则似然为
MRFx的联合分布可表示为Gibbs函数[9],
式 (8) 中, Z是规范化常量, Vc (x) 表示派系c的潜力, C是根据邻域系统确定的所有派系的集合, T是温度参数。本文使用Potts模型表示派系潜力, 运用式 (7) 和式 (8) 到式 (4) , 并对其进行负对数变换, 得到
根据模拟退火MCMC方法, 为温度参数T定义一个冷却进度表,
式 (10) 中, i=1, 2, …, I表示MCMC算法迭代的数目, C是冷却因子, 本文设置T (0) =4, C=0.97。若给定一幅脑部MR图像y和标签x (0) 的初始配置, 则可计算出模型参数。定义用来表示从x (i) 随机移动的跳跃密度Q (.|x (i) ) 符合高斯分布, 每次迭代中, 从建议密度Q (x* (i+1) |x (i) ) 提取一个候选x* (i+1) , 从均匀分布u (0, 1) 提取一个随机序列, 计算每个体素j的接受率
如果uj<αj, 接受模拟xj (i+1) =xj* (i+1) , 否则拒绝它并保持类标签与上一次迭代xj (i+1) =xj (i) 相同, 当达到最大迭代次数时停止, 如算法1所示。
算法I:体素分类的MCMC采样
2.2 偏场校正
在利用MCMC体素分类之后, 可以获得分割结果x*和最小能量E={Ekj, xj=k, j∈S}。归一化后验概率n={nkj;j=1, 2, …, N, k=1, 2, …, K}作为软分割结果
根据理想MR图像的分段常数性质, 定义软分割与对应平均μk的积作为存储的图像
本文使用奇异值分解 (SVD) 求解下列最小二乘拟合问题, 以此估计偏场。
式 (14) 中, ./表示点对点划分, 根据估计的组合系数, 获得偏场
偏场损坏的图像可恢复如下:
2.3 CSA进行参数估计
第三步是通过最大化后验概率p (Θ|y (t) , x (t) ) 来学习当前图像强度y (t) 和配置x (t) 给出的最优参数
式 (17) 中, p (Θ) 是参数的先验概率, 这个先验指的是基于马尔可夫性质信息p (Θkj) =p (xj=k|kj) 的体素, 可以通过MRF能量计算得到, 定义它作为这些项的混合来平衡参数的收敛和多样性, 对于每个体素j∈S
式 (18) 中, v是平衡常量, Πkj=πk表示分体素全球先验。给定任意具体参数集Θ, 即可计算式 (18) 中所示优化问题的目标函数。为了实现全局最优, 采用CSA[10]求解该问题, 以群体方式模拟所有可能参数。CSA是一种进化优化算法, 通过迭代生成一群编码抗体来寻找全局最优解。本文中, 抗体群np设为100, 定义每个抗体为一个候选参数集Θ, 将有特定抗原的抗体Θk的亲和力定义为后验似然p (Θ|y, x) , 迭代优化过程由下列六个主要步骤组成:
第一步:评估每个抗体的亲和力, 根据其亲和力按降序排列所有抗体;
第二步:从当前群体中选择Ns个抗体, 克隆它们形成克隆群。对于亲和力排序为j的抗体, 定义其克隆的数目正比于其亲和力排序, 如式 (19) 所示。
式 (19) 中, β是常量, 用来控制克隆率, round (.) 用来将实数变换到与其最接近的整数;
第三步:分别对概率为phm和pre的克隆群运用超突变和受体编辑操作。超突变是在动态范围±10%内随机改变抗体的值, 目的是局部搜索最优解。受体编辑是在动态范围±100%内随机改变抗体, 实现全局搜索;
第四步:评估克隆群中抗体的亲和力, 根据其亲和力按降序排列;
第五步:选择克隆群中排名靠前的抗体代替记忆细胞集中较低亲和力的40%抗体, 保证记忆细胞集保存迄今为止获得的最优解, 以便最高亲和力的抗体按代递增;
第六步:用随机生成的抗体代替剩余集中具有最低亲和力的10%抗体, 对新群体引入多样性。
重复迭代这个过程直到达到最大迭代数目, 如图1所示。
2.4 总结
给定K-平均算法产生的初始分割结果, HMRF-CSA算法迭代执行基于MCMC的体素分类、偏场校正和基于CSA的模型参数估计, 直到算法收敛。一旦达到收敛, 则获得最终分割结果、偏场和模型参数。HMRF-CSA算法的主要步骤见算法II。
算法II:HMRF-CSA脑部图像分割算法
3 实验结果
本文从Brain Web数据集[11]获取仿真TI加权脑部MR图像, 比较提出的HMRF-CSA算法与现有的e HMRF算法、GAMIXTURE包中GA-EM算法、D-C算法、FSL包中的经典HMRF-EM算法和SPM包中的统一分割程序。Brain Web数据集提供的一组仿真脑部图像, 这些图像具有各种INU和噪声级别的解剖模型仿真, 每个仿真研究的维度为181×217×181, 体素大小为1 mm×1 mm×1 mm。
图2分别显示了仿真图像中具备40%INU和7%噪声的第88个横切片, 偏长矫正图像, 估计的偏场, 使用六种算法获得的分割结果和地面实况组织图。其中图2 (a) 表示仿真图像的第88个横切片 (7%噪声和40%INU) ; (b) 表示INU校正图像; (c) 估计的INU; (d) 表示HMRF-EM算法的结果; (e) 表示D-C算法的结果; (f) 表示SPM算法的结果; (g) 表示GA-EM算法的结果; (h) 表示e HMRF算法的结果; (i) 表示HMRF-CSA算法的结果; (j) 表示地面实况。可以看出, 本文算法产生的分割结果比其他算法产生的结果更接近地面实况。
接下来, 在两组仿真MR图像上对这些算法进行进一步比较。第一组MR图像包含有20%INU和噪声级别范围从1%到7%的四个图像, 使用骰子相似度系数 (DSC) [12]定量评估每个脑部组织类型分类的性能。
式 (20) 中, Vs (k) 是分割结果中脑部组织类k的体, Vg (k) 是在地面实况上对应的体, |V|代表体V中体素的数目。通过正确分类的脑部体素百分比来计算分割精度, 并评估整体精度。图3表示六种算法获得的分割精度。
从图3可以看出, 在大部分仿真图像中, 本文算法在划分每个脑组织和分类整个脑部体方面都具有较高的精度。而且, 随着噪声和INU级别的增加, 提出算法的精度下降幅度比其他算法低, 这表明本文算法具有较强的抵制噪声和INU影响的能力。
第二个测试组包含有40%INU和噪声级别范围从1%到7%的四个图像, 六种算法获得的分割精度如图4所示。
从图4可以看出, 本文算法能在高噪声和INU级别下保持良好分割性能。
4 讨论
4.1 参数设置
本文提出的HMRF-CSA算法中, 有三组需要近似的参数, 包括MCMC推断、INU估计和基于CSA的参数近似。在INU近似中, 以正交多项式的阶来权衡考虑近似精度和计算复杂度, 由于INU变化非常慢, 对于INU近似, 10个三阶多项式已经足够。式 (19) 中权重参数v决定MRF先验, 较大的v能使MRF的作用更大, 另一方面, 小v则更支持GMM先验。CSA本身需要很多参数, 文献[13]对此进行了详细讨论。本文使用CSA程序的经验参数设置:群体大小Np=100、记忆集大小Nm=0.3Np、选定抗体的数目Ns=0.5Np、克隆率常量β=0.5, 超突变概率phm=0.8、受体编辑概率pre=0.1和最大的代Nt=20。
4.2 计算复杂度
计算机程序的性能与许多因素有关, 包括计算机处理能力、数据表示、编程语言和编码实现等[14]。本文评估了HMRF-CSA算法的计算复杂度, 本文算法在每次迭代中顺序执行MCMC推断、偏场估计和基于CSA的参数估计。设定对于有N个体素的一副图像, MCMC推断的计算复杂度为O (N) ;偏场估计仅进行一些矩阵计算;基于CSA的参数估计的复杂度为O (Np+NcK) , 其中Np表示群体大小, Nc表示总克隆数目。提出的迭代分割算法的迭代次数达到wmax后停止, 其线性整体计算复杂度O (N+Np+NcK) 。需要注意的是, MCMC方法的主要缺点是需要大量仿真图[15], 然而, 由于CSA为MCMC方法配置了一个良好的开始状态, 所以本文算法不需要许多仿真图。同时, MCMC方法的输出使CSA在有限代数之后成熟, 因此, 尽管本文HMRF-CSA算法涉及耗时的MCMC和CSA程序, 然而, 其计算复杂度只稍微高于传统分割方法。
5 结论
提出了HMRF-CSA脑部MR图像分割算法, 在基于HMRF模型分割中结合了CSA和MCMC, 本文算法能够有效的用于基于HMRF模型估计的图像分割问题。在仿真脑部MR图像上进行实验, 将本文算法与GA-EM算法、D-C算法、SPM和FSL软件包算法进行比较, 实验表明该算法获得了更好的分割精度。
未来会将本文算法应用于其他的脑部MR数据库上, 并与其他新颖技术相结合, 进行大量的实验, 继续提高分割精度。
马尔可夫 篇2
In this paper,we define a model of random dynamical systems(RDS)on graphs and prove that they are actually homogeneous discrete-time Markov chains.Moreover,a necessary and sufficient condition is obtained for that two state vectors can communicate with each other in a random dynamical system(RDS).
作 者：郑洁 刘朝阳 ZHENG Jie LIU Chao-yang 作者单位：郑洁,ZHENG Jie(Department of Applied Mathematics,Donghua University,Shanghai 20,China)
刘朝阳,LIU Chao-yang(Zhengzhou Dongxifang Computer Network Engineering Limited,Henan University,Zhengzhou 450008,China)
马尔可夫信源极限熵求解方法解析 篇3
1 马尔可夫信源
⑴马尔可夫信源。我们说实际的信源一般都是有记忆的信源, 而且这种有记忆信源在任一时刻发出符号的概率通常只与前面若干个符号有关, 而与更前面的符号无关, 因此我们可以认为信源在某一时刻发出的符号与信源的状态有关。若信源输出的符号序列和状态序列满足下述的两个条件:某一时刻信源的输出仅与信源的当前状态有关;信源的状态只由当前的输出符号和前一时刻信源状态唯一确定。我们称这样的信源为马尔可夫信源。
⑵马尔可夫信源的极限熵。若信源以长度为N输出符号序列, 则信源的平均符号熵为其中HN (X) 是信源的矢量熵。当N→∞时, 此时称为信源的极限熵, 极限熵是真正描述实际信源熵的表达方式。它规定了平稳离散有记忆信源输出符号序列中平均每个信源符号的熵值, 代表了一般离散有记忆信源平均每发出一个符号所提供的信息量。事实上, 当信源记忆长度很长, 趋于无穷大的时候, 要计算联合熵或极限熵是很困难的, 它需要测定信源的无穷阶联合概率和条件概率, 这是很难达到的, 因此, 我们在实际计算时, 我们往往只考虑有限记忆信源的熵, 用有限的条件熵或平均符号熵作为极限熵的近似值。
由此可见, 当信源是有记忆m阶的马尔可夫信源时, 我们用条件熵作为极限熵的近似值。而求解条件熵的关键就是要得到马尔可夫信源稳定后 (N→∞) 各个状态的极限概率。
2 马尔可夫信源极限熵求解案例分析
极限熵并不是在任何情况下都存在。通常, 对于一个n元m阶的马尔可夫信源, 只有在平稳状态下, 各个状态的状态极限概率都存在时, 才可以计算出极限熵。因此, 求解马尔可夫信源的极限熵关键在于如何求解马尔可夫信源稳定后各个状态的极限概率。下面, 就以一个案例来说明马尔可夫信源的极限熵的求解方法。
举例:设有二元2阶马尔可夫信源, 其原始信源X的符号集为{x1=0, x2=1}, 其状态转移图如下图所示, 求该马尔可夫信源的极限熵。
解: (1) 首先求解各状态极限概率
这样二元2阶信源X∈{0, 1}得到的状态空间{e1, e2, e3, e4}和相应的一步转移概率构成的2阶马尔可夫信源模型为:
将一步转移概率代入上式得方程组:
求求解解以以上上方方程程组组, , 则则可可算算出出该该信信源源的的状状态态极极限限概概率率为为:
(2) 求解该信源的极限熵
计算出该马尔可夫信源的状态极限概率后, 根据状态转移图提供的一步转移概率, 我们就可以计算出这个2阶马尔可夫信源的极限熵了。
3 总结
马尔可夫信源属于有记忆信源, 当信源在某个状态时, 只取决于这个时刻发出的符号与之此时刻之前的符号状态有关。马尔可夫信源不同于一般有记忆信源之处在于它用符号之间的转移概率来描述符号间的关联性, 即马尔可夫是以转移概率发出每个信源符号的。计算马尔可夫信源的极限熵时, 首先要考虑该信源的极限熵是否存在, 若极限熵存在, 则需先计算出该信源各个符号状态的极限概率, 再根据极限概率和转移概率求出极限熵。通过计算极限熵, 我们可以计算出信源存在的冗余度, 为信源的编码奠定基础。
参考文献
[1]陈运.信息论与编码.北京:电子工业出版社, 2009.
[2]张旭东, 等.图像编码基础和小波压缩技术.北京:清华大学出版社, 2004.
[3]吴家安, 等.语音编码技术及应用.北京:机械工业出版社, 2006.
[4]钟家恺.通信原理教程.北京:科学出版社, 2003.
马尔可夫 篇4
光线追踪算法用于计算机生成真实感图像, 在电影, 电视和辅助设计领域有着广泛的应用。光线追踪算法和另一大图像生成算法光栅化相比, 其图像生成的质量和真实性均占优, 但是速度要比光栅化算法慢很多。自从该算法被提出以来, 研究人员提出了多种对此算法的改进, 比如光子映射算法和辐射度算法。这些算法和传统光线追踪算法比起来在生成图像质量相同的情况下速度快了很多。但是这些算法只针对特定的场景布置才能起到加速效果, 在一些复杂场景特别是有大量光滑和光泽表面的情况下会出现大量误差, 并且计算速度会大幅下降。
光线追踪算法的本质是求解光传输方程。通过光能守恒定律我们可以把积分空间离散的传输方程转换为场景中面积度量的积分。求解这个积分的方法是把它转换为一个无穷维的高阶积分, 并用蒙特卡罗方法采样整个积分空间, 使积分结果按概率趋近于方程解。由于渲染场景的复杂性和光照模型的多样性, 求解这个问题的最大难点在于如何在有限的时间内使蒙特卡罗方法逼近真实解。目前, 光线追踪主要还是应用于离线渲染领域。如何在保证图像质量的情况下提高算法速度就成了亟待解决的问题。
二、问题描述
2.1光传输方程。
真实感渲染中, 核心算法之一是模拟光从光源发出, 在场景物体和介质之间进行交互, 传播, 最终被吸收或者进入摄像机的过程。从物理学的角度来看, 这个过程就是电磁波 (辐射) 与实物粒子的碰撞与反射。我们将一个表面上一点p处ω0方向上辐射出射能量定义为L0 (p, ω0) , 那么根据热辐射定律, 出射能量可以表示为入射能量在表面BRDF函数和方向ω的积分, 该积分可以用下列公式表示:
其中Le (p, w0) 是表面发出的辐射;f (p, ω0, wi) 是表面的BRDF函数, 表征物体表面对辐射的反射和吸收特性。积分空间s2是场景中的所有表面。由于积分中含有即立体角, 为了能计算这个积分, 我们希望它能统一到积分空间中。有限空间的场景中, 沿着射入一个表面的光线能找到另一个场景中的点。于是我们定义函数, 它表示点p处沿着方向找到的另一个场景中的点。如果这样的点不存在, 则t函数返回一个特殊的值。将函数t带入式 (1) 可以得到
(2) 式与 (1) 式相比, 右侧的积分函数和左侧统一, 但是依然无法直接进行计算, 因为它依然含有ω。考虑到我们在上一步的变换同样适用于L (p, ω0) , 为此我们定义
其中P'为沿着点p的出射光线遇到的场景中的点。由于场景中存在遮挡关系, 我们定义几何函数
函数表示点p和P'的可见性。如果两点之间没有遮挡, 那么V的值是1, 否则是0.将上述函数引入式 (2) 中, 我们得到了光传输方程的表面形式:
注意到上式中左边和右边都有函数L。现在我们考虑一条从光源到摄像机的完整光路径, 我们把P0处的式 (3) 带入P1处的式 (3) 得到
如果光路径上还有其他点, 也可以同样按照上述方法带入, 从而把光传输方程转化成路径形式。有了式 (4) 就可以在某些简单场景中计算数值解了。但是对于大多数复杂场景, 光路径的长度可能很长甚至是无限长, 式 (4) 就变成了一个无穷维积分, 这样的积分是无法计算数值解的。对于这样的场景, 我们用2.2中的采样算法进行计算。
2.2 Metropolis-Hastings采样算法。
其中表示以概率a接受从X状态到X’状态的变化, 否则不变。由 (5) 式得
以这种采样的方式计算积分的话, 就是MCMC算法:
其中g (x) 可以是任意的函数, xi服从分布.
式 (4) 可以简化写为
使用MCMC算法, 式 (6) 可以用下式计算:
三、算法流程与实现
本文算法的实现使用了开源项目mitusba render作为输入和输出, 核心算法流程图如下:
四、实验与结果分析
实验使用光线追踪测试场景veach, 生成的图像如下:
通过生成的图像我们可以看到, 使用MCMC算法生成的图像整体具有较好的真实感, 并且能够表现出光线多次反射对场景中物体的影响。但是, 该算法的速度比较慢, 在2450M处理器上渲染上图使用了15分钟26秒。如何提高算法的性能将是一部研究的重点。
参考文献
[1]Veach.Robust Monte Carlo methods for light transport simulation[D].斯坦福:斯坦福大学, 1997.
[2]Jakob W., S.Marschner.Manifold exploration:A Markov Chain Monte Carlo technique for rendering scenes with difficult specular transport[J].ACM Trans, Graph, 2012, 31, 4, 58:1-58:13.
[3]Jaakko L., K.Tero, L.Samuli, A.Miika, D.Fr'edo, A.Timo.Gradient-Domain Metropolis light transport[J].ACM Trans.Graph, 2013, 32, 4, 95:1-95:11.
[4]Jakob W.2012 Mitsuba v0.4 EB/OL].http://mitsuba-renderer.org.
[5]Jaakko L., K.Tero, L.Samuli, A.Miika, D.Fr'edo, A.Timo.GradientDomain Metropolis light transport[J].ACM Trans.Graph, 2013, 32, 4, 95:1-95:11.
[6]Jakob, W.Light transport on path-space manifolds[D].伊萨卡:康奈尔大学.2013.
马尔可夫 篇5
关键词：隐马尔可夫模型,学习行为,滑动窗口,抄袭,学习异常
0 引言
网络教学打破传统课堂教学的时空限制, 具有便捷、泛在访问等特点, 正成为时下高校流行的教学方式。学生可通过丰富的网络媒介 (如浏览网页、讨论板、博客、wiki、提问等) 进行学习。当前已有学者使用BP神经网络和C4.5算法对学生使用网络媒介的学习行为进行效果评估[1,2]。为了强化学生对知识的灵活运用, 教师通常会引入虚拟实验实训、在线测试之类的系统对学生进行知识训练。文献[3]使用隐马尔可夫模型HMM (Hidden Markov Model) 实现知识点智能引导。笔者曾使用模糊聚类方法对CSCL学习者的混合分组进行了基础研究[4]。在小组学习过程中, 有时候学生为了应付教师, 往往弄虚作假。对于学生抄袭或学习异常的学习行为, 系统反馈的有可能是虚假信息, 教师难以及时发现隐藏在虚假信息背后的真实行为, 无法及时帮助学生纠正错误。
隐马尔可夫模型在语音识别、网络异常检测、行为异常检测方面得到广泛的应用[5,6]。为了从大量学习记录中挖掘和提取有用信息, 本文引入隐马尔可夫模型, 以数据库课程SQL在线测试系统为例, 建立基于隐马尔可夫模型的学习行为评估模型, 评估和发现学生的抄袭和异常行为。
1 隐马尔可夫模型
1.1 隐马尔可夫模型定义
隐马尔可夫模型是一个双重随机过程, 包含两个随机变量序列:一个是观察不到的马尔可夫链, 用来描述状态的转移, 用转移概率表示;另一个是可以观察到的随机序列, 用来描述状态与观察值的关系, 用观察值概率表示。完整的隐马尔可夫模型由一个五元组λ= (S, V, A, B, π) 表示[7], 其中:
(1) S为隐藏状态集合, S={s1, s2, …, sN}, |S|=N, 并记t时刻的状态为qt, qt∈S。
(2) V为观察符号集合, V={v1, v2, …, vM}, |V|=M并记t时刻观察到的符号为ot, ot∈V。
(3) A为状态转移概率矩阵, A= (aij) , aij表示如果在t-1时刻状态为si, 则在t+1时刻转移到状态sj的概率, 即aij=P (qt+1=sj|qt=si) 1≤i, j≤N。
(4) B为状态的观察符号概率分布, B={bj (k) }, bj (k) 表示在状态sj下观察到符号vk的概率, 即bj (k) =P (ot=vk|qt=sj) 1≤k≤M, 1≤j≤N。
(5) π为初始状态的概率分布, π={πi}, πi表示在时刻t=1时, 处于状态si的概率, 即πi=P (q1=si) 1≤i≤N。
隐马尔可夫模型的性质完全由A、B、π所确定, 为了方便, 简记为λ= (A, B, π) 。如无特别说明, 这里指的是一阶隐马尔可夫模型。
1.2 隐马尔可夫模型的应用
已知观察序列O= (o1, o2, …, oT) 和模型λ= (A, B, π) , HMM在实际的应用中要解决下列三个问题:
(1) 评估问题:求模型λ产生观察序列O的条件概率P (O|λ) , 可使用前向算法求解。
(2) 解码问题:求模型λ产生观察序列O的最可能的状态序列, 可使用Viterbi算法求解。
(3) 学习问题:使用观察序列O, 调整模型λ参数, 使得条件概率P (O|λ) 最大, 可使用Baum-Welch算法求解。
1.3 Baum-Welch算法介绍[8]
给定模型λ和观察序列O:
(1) 定义前向变量:
表示在给定模型λ的条件下, 在t时刻, 产生部分观察序列 (o1, o2, …, ot) , 并处于状态si的概率。前向变量可由下式进行迭代计算:
其中, α1 (i) =πibi (o1) 。
(2) 定义后向变量:
表示t时刻状态为si的条件下, 在t+1时刻到最后, 产生部分观察序列 (ot+1, ot+2, …, oT) 的概率。后向变量可由下式进行迭代计算:
其中, βT (i) =1。
(3) 定义输出概率:
(4) 定义t时刻在状态si, t+1时刻在状态sj的概率:
(5) 定义t时刻处于状态si的概率:
(6) 定义从状态si转移到状态sj的期望次数与从状态si转移的期望次数之比:
(7) 定义处于状态sj且观察符号为vk的期望次数与处于状态sj的期望次数之比:
(8) 定义t=1时刻处于状态si的概率:
构成新的HMM重估模型。给定初始模型λ= (A, B, π) , 利用训练序列, 通过反复迭代, 计算重估模型。当重估模型收敛时, 即可得出局部最优的HMM重估模型参数。Baum-Welch算法如下:
输入:A, B, π, O
输出:
(1) 初始化终止条件δ。
(2) 迭代计算:
(3) 如果Δp>δ, 继续迭代;否则终止。
2 建立隐马尔可夫模型
2.1 SQL在线测试系统介绍
本文设计了一个SQL在线测试系统。学生登录系统, 在题库中选择训练题, 按要求编写并在线提交SQL语句。测试系统编译并运行SQL语句, 并对运行结果进行自动评估。测试系统反馈六种结果:结果正确、SQL错误、行数不等、列数不等、数值不等和列名不等。学生做题是随机的, 每道题都会产生6种反馈结果中的一种。教师无法直接了解学生对题目知识的掌握程度, 只能通过学生做题时系统的反馈结果进行观察。
笔者把SQL在线测试系统应用在计算机专业二个教学班上, 经过一学期的实际教学使用, 共收集到94个学生的2万多个做题记录。
2.2 设置模型参数
学生选题是随意的、无序的, 但做题尝试是和时间相关的。学生的每次做题尝试, 都可看作是对知识点掌握程度的一个反映 (观察值) 。做题尝试序列反映知识点掌握程度状态的转移。
根据学生的学习情况, 定义对知识点的掌握程度状态集S={完全掌握, 基本掌握, 了解一点, 完全不会, 抄袭}, 状态数N=5。定义观察符号集V={结果正确, SQL错误, 行数不等, 列数不等, 数值不等, 列名不等}, 观察符号数M=6。初始为全体学生做题序列里符号vi的数量。根据学生的学习状态, 估计初始状态转移矩阵和状态观察符号概率分布, 作为模型的初始训练参数。
2.3 模型训练
以学生为单位, 以时间先后次序, 生成每个学生的做题反馈结果序列 (观察序列) 。先从全体学生样本中随机抽取15%样本作为未知集 (待检验集) 。然后根据教师观察学生的课堂练习情况, 通过人工方式把剩下的学生样本归为:正常集、抄袭集和异常集。正常集中抽取65%数据作为训练集, 剩下35%作为基准集。
本文前面列出的Baum-Welch算法公式只适用于单个训练样本, 对于训练集里存在多个学生样本序列的情况, 这里使用适用于多训练样本的推广的Baum-Welch算法。对训练集应用推广的Baum-Welch算法进行训练, 得出正常行为HMM模型。
3 模型应用
3.1 引入滑动窗口技术
由于每个学生的做题记录是不同的, 因而观察序列的长度也差别较大。因此不能简单地直接比较各个学生观察序列的模型输出概率。为了使输出概率的比较更具科学性, 引入滑动窗口技术[9]来分割观察序列。令l为滑动窗口大小, 滑动窗口每次向后移动一位, 整个观察序列可分为T-l+1个子序列 (T为整个观察序列的长度) 。令W={wk} (1≤k≤T-l+1) 为子序列的集合, 使用滑动窗口求解所有子序列输出概率P (wk|λ) 的滑窗-前向算法如下 (为了突出输出概率, 方便比较, 这里取输出概率的对数) :
输入:N, O, T, π, A, B, l
输出:log P (wk|λ)
(1) 循环:wk= (ok, ok+1, …, ok+l-1) k=1, …, T-l+1
(2) 对每个wk进行:
(3) 初始化:α1 (i) =πibi (ok) 1≤i≤N;
(4) 迭代计算:
(5) 迭代终止:
(6) 循环结束。
3.2 评估学生做题序列
使用滑窗-前向算法应用正常行为HMM模型对基准集每一个样本求所有子序列的输出概率对数均值:
取所有样本的均值作为学生做题的正常值。
对抄袭集和异常集, 分别对集合里每个学生样本的观察序列使用滑窗-前向算法求各个子序列的log P (wk|λ) 。将所有的子序列标记为“抄袭”。定义所有“抄袭”的子序列数φhigh与总子序列数φ的比值为抄袭度:
将所有的子序列标记为“异常”。定义所有“异常”的子序列数φlow与总子序列数φ的比值为异常度:
4 实验结果
4.1 实验数据与分析
由于学生做题的平均尝试次数约为7次, 这里设定滑动窗口大小为7。正常集、抄袭集、异常集里所有样本的对比如表1所示。
由表1可知, 抄袭集的输出概率均值比基准集高, 原因是学生都是抄袭别人正确答案, 观测值都是“结果正确”居多, 其它观测值很少, 自然输出概率均值比正常值高;异常集的输出概率均值比基准集低, 原因是学生正确做题的数量不多, 但是错误反馈比基准集多。
对各个数据集的每个学生样本计算其抄袭度和异常度, 结果如表2所示 (限于篇幅, 只列出部分典型样本数据) 。
由定义可知, 正常的抄袭度和异常度应该都是在0.5上下小幅浮动。表2的数据显示, 抄袭集的抄袭度和异常集的异常度都高于正常值。通过人工分类的抄袭集和异常集, 在HMM里得到正确标记。
HMM对未知集样本评估结果显示:S128比值偏高, 属典型的抄袭行为;S137和S244显示正常;S235比值偏低, 可能学习存在问题。通过课堂观察和学生面谈等人工调查, 结果显示, S128课堂常开小差, 大部分题目抄袭别人答案;S137独立完成;S235对知识掌握得不好, 常常出现语法等错误;S244是学习尖子, 喜欢独立钻研。观察S244的学习序列, 发现某些错误连续多次出现, 长度是普通学生的2倍多, HMM因而把其判定异常。除S244估计有偏差外, HMM对未知数据集其余样本的学习行为评估基本正确。
4.2 HMM和其它方法对比
当前还没发现和本文相似的研究文献, 为了对比HMM对各个数据集评估的正确率, 这里使用一种简单可行的方法:以学生每道题的平均尝试次数 (题均数) 来检测抄袭和异常学习行为。正常的题均数的范围为3~9之间, 低于此区间的定为抄袭, 高于此区间的定为异常。HMM与题均数法的结果对比如图1所示。
由图1可以看出, 通过人工分类的抄袭集和异常集, HMM都能正确标记, 但题均数方法就存在偏差, 对未知集的评估比HMM差。总的来说, HMM方法的正确率比题均数法要高。
5 结语
本文介绍了隐马尔科夫模型的原理。以SQL在线测试系统为例, 使用HMM建立学习行为评估模型。引入滑动窗口技术来解决学习序列长度不一而影响输出概率的问题。使用正常的学习序列进行训练, 由此评估学生的学习行为是否偏离正常。实验结果表明, HMM模型识别抄袭和异常学习行为的正确率比普通方法要高。虽然以SQL在线测试系统为例, 但本文使用的方法完全可以推广到其它一般系统。本文在实际应用中也存在一些需要改进的地方, 如对知识点掌握程度的状态集的分类有待更深入研究, 可使用半监督的建模方法[10]改进人工标记数据集进行建模。
参考文献
[1]姜华, 赵洁.基于BP神经网络的学习行为评价模型及实现[J].计算机应用与软件, 2005, 22 (8) :89-91.
[2]范洁, 杨岳湘, 温璞.C4.5算法在在线学习行为评估系统中的应用[J].计算机工程与设计, 2006, 27 (6) :946-948.
[3]翟琳琳, 陈仪香.隐马尔科夫模型在智能学习系统中的应用[J].计算机工程与应用, 2007, 43 (6) :178-180.
[4]黄志成.基于模糊聚类的CSCL学习者混合属性分组[J].计算机应用与软件, 2011, 28 (2) :118-121.
[5]温凯, 郭帆, 余敏.自适应的Web攻击异常检测方法[J].计算机应用, 2012, 32 (7) :2003-2006.
[6]李战明, 宋丙菊.基于隐马尔可夫模型的ATM机用户异常行为识别[J].兰州理工大学学报, 2012, 38 (5) :77-81.
[7]Rabiner L R.A tutorial on hidden Markov models and selected applications in speech recognition[J].Proceedings of the IEEE, 1989, 77 (2) :257-289.
[8]杜世平.隐马尔可夫模型的原理及其应用[D].四川:四川大学, 2004:1-10.
[9]张响亮, 王伟, 管晓宏.基于隐马尔可夫模型的程序行为异常检测[J].西安交通大学学报, 2005, 39 (10) :1056-1059.
马尔可夫 篇6
1 马尔可夫预测模型的建立
1.1 马尔可夫链的定义及性质
马尔可夫链数学表述:定义1:设有离散的随机变量过程{Xn,n∈N,N=0,1,2,L},X所有可能取值的全体称为{Xn}的状态空间,记为E={x1,x2,K}。若对任意正整数n及任意xi,(xi1xi2,L,xin,xin+1)∈E,只要:
P(X1=xi1,X2=xi2,L,Xn=Xin)>0 (1)
有:P(Xn+1=xin+1|X1=xi1,X2=xi2,L,Xn=xin)
=P(Xn+1=xin+1|Xn=xin) (2)
则称{Xn}为马尔可夫链[1,2,3,4]。
定义2:若{Xn}为马尔可夫链,对任意xi,xj∈E;m≠n,总有:
P(Xn+1=xj|Xn=xi)=P(Xm+1=xj|Xm=xi) (3)
则称{Xn}为齐次的马尔可夫链。对齐次的马尔可夫链,有对任意m,k∈T,有:
Pij(m;k)=Pij(k),i,j∈E (4)
式中,Pij(m;k)为系统在m时刻处在状态i,经k步状态转移到达状态j的概率;Pij(k)为系统从状态i,经k步状态转移到达状态j的概率,此时转移概率与初始时刻无关,k取1时,Pij记为Pij(1)。齐次马尔可夫次链完全由其初始分布{P0(i0),i0∈E}及其一步状态转移概率矩阵p=(Pij),i,j∈E所决定。若已知时刻n的绝对分布p(n)={pn(j),j∈E},则时刻n +1的绝对分布为:
undefined
1.2 模型的建立过程
假定某一种被预测的事件有E1,E2,…,En,共n个可能的状态。记Pij为从状态Ei转为状态Ej的状态转移概率,作矩阵:
undefined
则称P为状态转移概率矩阵。为了运用马尔可夫预测法对事件发展过程中状态出现的概率进行预测,还需要再介绍一名词:状概率πj(k)。πj(k)表示事件在初始(k=0)时状态为已知的条件下,经过k次状态转移后,第k个时刻(时期)处于状态Ej的概率。根据概率的性质,显然有:
undefined
从初始状态开始,经过k次状态转移后到达状态Ej这一状态转移过程,可看作是首先经过(k-1)次状态转移后到达状态Ei(i=1,2,…,n),然后再由Ei经过一次状态转移到达状态Ej。根据切普曼—柯尔莫哥洛夫方程有:
undefined
若记行向量π(k)=[π1(k),π2(k),L,πn(k)],则由式(8)可得逐次计算状态概率的递推公式:
undefined
式中,π(0)=[π1(0),π2(0),L,πn(0)]为初始状态概率向量。
2 文登市降水量预测
2.1 研究区概况
文登市位于山东半岛东部,36°52′—37°23′N、121°43′—122°19′E;西阻于昆嵛山,与烟台市牟平区和乳山市相接,北连威海市环翠区,东邻荣成市,南濒黄海。文登市总面积1645km2,海岸线155.88km;地处北温带,属于大陆性季风气候,四季分明;年均气温11.5℃,有气象资料记录以来的极端最高气温为36.4℃,极端最低气温-25.5℃;降水量762.2mm,降水分布不均,夏季较为集中,6—9月份降水量约占全年70%,春秋季降水偏少,常发生干旱;年日照时数2540.7h,无霜期194d;冬季漫长,盛行从大陆北部吹来的干冷冬季风,气温偏低,为半岛地区低温站点;夏季最短,盛行从海洋吹来的暖湿夏季风,春秋两季属冬夏季风转换期,春季大风天气较多。
2.2 模型的应用
以山东省文登市降水量作为研究实例,文登市1953—2005降雨量的变化情况见图1。考虑到山东省文登市年降水量变化的3个状态,即“枯水年(≤550mm)"、“平水年"和“丰水年(>800mm)",分别记E1为“枯水年"、E2为“平水年"、E3为“丰水年"。文登市1953—2005年年降水量的状态变化情况见表1。
由表1计算可得该地区年降水量的状态转移概率矩阵为:
undefined
根据建模过程可知,如果某一事件在第0个时刻(或时期)的初始状态已知,即π(0)已知,则利用递推公式(9)就可求得它经过k次状态转移后在第k个时刻(时期)处于各种可能状态的概率,即π(k),从而得到该事件在第k个时刻(时期)的状态概率预测。如果将2005年的年降水量状态记为π(0)=[0,0,1](2005年处于“丰水年"状态),则将状态转移概率矩阵(10)式及π(0)代入递推公式(9)式,就可求得2006—2010年可能出现的各种状态的概率,见表2。
2.3 城市雨水的利用
由于文登市属于大陆性季风气候,而且降水夏季较为集中,春秋季降水偏少,因此有必要利用城市雨水资源解决春秋季的水资源不足,实现水资源利用的可持续性发展。由图2和表2可知,未来的5年来是丰水年的概率很高,因此协调春秋季与夏季之间水量不平衡则极其重要。实施城市雨水收集工作比农田集水的形式便捷,收集利用率高,一般可达到75%以上。这主要是城市的大多数街道为柏油路面,雨水通过改造后的地下管网汇集起来,经过处理后再使用。随着城市化发展,城市生活用水、工业用水和生态用水大幅度增加,不仅挤占农业用水,也加重了地下水的负担。在长期的发展规划中,如何合理、可持续地利用水资源具有战略性意义,因此城市雨水资源利用是解决城市水危机的必然选择。
许多地区对城市雨水利用未引起足够的重视,其原因一是观念落后,几千年以来人类对雨水放任自流已形成根深蒂固的思维定势;二是经济账不清,对城市水危机和雨水资源化认识不足,未核算雨水利用的投入产出;三是缺乏政府行为,规划、立法、规范等问题亟待解决,因此建议:①加强宣传,确立雨水是宝贵的水资源的共识。通过各种媒体、社区活动等形式进行宣传,提高全社会的雨水资源利用意识,取得民众的广泛参与和支持,不能让雨水资源放任自流地白白浪费,要确定综合利用在前、排放在后的指导思想。②结合规划,改建雨水专用排水系统。现在大多数城市排水管网都是雨水和污水合流,雨水比污水的水质好,易处理,因此将雨水和污水分流处理是较好利用雨水资源的一个前提条件。传统的排水系统主要是将雨水及时排走,避免路面积水,而雨水利用则是需要将雨水进行存储、运输,所以需要把现有的排水系统改造成既能很好利用雨水又不影响城市环境的雨水专用排水系统。③制定法规,鼓励和强化对雨水的利用。针对研究区的水资源和环境状况,并借鉴发达国家的经验,充分合理地利用雨水资源,应尽早制定一系列有关雨水利用的法律法规和条例,规范新建城区、居民小区和工、商业区,设计雨水利用及蓄水设施等内容。同时,采用鼓励对雨水的利用政策,调动开发商和企事业单位的积极性。④把城市雨水资源的利用纳入城市整体规划。城市雨水资源的利用是一项造福子孙后代的系统工程,应纳入城市整体规划。因此,在统一考虑城市雨水利用与城市建设、水资源优化配置、生态建设时,应把集水、蓄水、处理、回用、入渗地下、排水等纳入城市建设规划之中。
3 结论
在水文及气象科学中,降水量的预测是一项非常重要的工作。本文所提出的预测方法具有如下特点:①由表2和图2可知,未来的5年山东省文登市的降雨量出现"枯水年"的概率最低,出现“平水年"的概率居中,出现“丰水年"的概率最高。总之,未来5年山东省文登市是丰水年的可能行较大。②预测结果为降水量的某一个状态的发生概率,而不是一个具体的数值,在可以完全满足实际工作的前提下,预测的范围更广,其可靠性也得到了提高。③根据最后计算出来的状态概率分布,求得所预测的降水量的具体值仍是一个有待解决的问题,马尔可夫链方法为这种预测提供了一个新的思路。
摘要：在介绍马尔可夫链的预测方法与模型的基础上,以山东省文登市为例,根据1953—2005年的降水量资料,对年降水量进行了状态分级;应用马尔可夫模型对该地区的降水情况进行了预测和分析,提出了对雨水资源利用的看法。结果表明,该方法客观、准确、可靠、简便,为区域降水的中短期预测提供了新的解决途径。
关键词：马尔可夫,降水量预测,雨水资源利用
参考文献
[1]王家生,刘嘉.随机过程基础[M].天津:天津大学出版社,2003.
[2]冯强,王昂生.用马尔可夫模型预测长江中下游地区旱涝灾害的试验[J].中国减灾,1997,7(4)∶17-20,24.
[3]丁晶,邓育人.随机水文学[M].成都:成都科技大学出版社,1988.
[4]张汉雄.用马尔可夫链模型预测宁南山区旱情[J].自然灾害学报,1994,3(1)∶47-54.
[5]宋志斌,张洪清,韩艳芳,等.对城市雨水资源利用的思考[J].水科学与工程技术,2006,22∶34-36.
[6]侯玉玲,张艳红,李春辉.城市雨水资源利用现状及发展建议[J].水科学与工程技术,2004,(6)∶11-13.
[7]郝红海,刘志强,刘畅.城市雨水资源利用研究[J].安徽农业科学,2007,35(1)∶163-164.
[8]夏乐天,朱元生生,沈永梅.加权马尔可夫链在降水状况预测中的应用[J].水利水电科技进展,2006,26(6)∶20-23,27.
马尔可夫 篇7
隐马尔可夫模型的概念是一般马尔可夫链概念的自然推广, 近年来, 在弱相依变量的建模上得到了广泛应用, 是研究发音过程、神经生理学与生物遗传等问题的有力工具。在理论方面, Leroox[1]与Bickel和Ratof[2]分别给出了隐马尔可夫模型在大数定律与中心极限定理方面的一些性质。在实际当中经常遇到隐藏链为非齐次马氏链的情况, 如动态的图形处理、气候的预测等均需要建立非齐次马尔可夫模型来处理[3,4]。所以研究隐非齐次马可夫模型的极限性质具有十分重要的意义[5,6]。
定义1.1 设S={1, 2, …, M}, T={1, 2, …, N}为两个有限集, {Xn, n≥0}与{Yn, n≥0}是概率空间 (Ω, F, P) 上的取值于S与T的随机变量序列.假设{Xn, n≥0}是非齐次马氏链其初始分布为 (q (1) , q (2) , …, q (M) ) , 转移矩阵为Pn= (an (i, j) ) M×M, i, j∈S, n≥1, 其中an (i, j) =P (Xn=j│Xn-1=i) , 称{Xn, n≥0}为状态链.它不能被直接观测到, 称为隐藏链;而能观测到的是{Yn, n≥0}, 称为观测链.如果存在矩阵B= (bij) M×N (i∈S, l∈T) 满足
P (X0=x0, Y0=y0, …, Xn=xn, Yn=yn) =q (x0) bx0y0a1 (x0, x1) …an (xn-1, xn) bxnyn。
则称{Xn, Yn, n≥0}为一个马尔可夫模型[7]。由于隐藏的马尔可夫链是非齐次的, 文中不妨称之为隐非齐次马尔可夫模型。
引理1[8] 设{Xn, Yn, n≥0}是如前定义的隐非齐次马尔可夫模型, f (x, y, z) 为定义在S×S×T上的实值函数, 令Fn=σ (Xm, Ym, 0≤m≤n) , 则有
引理2[9] 设{Xn, n∈N}是鞅差序列, {an, n∈N}为单调上升趋向于无穷的数列, 且
2 主要结果
定理1 设{Xn, Yn, n≥0}是如前定义的隐非齐次马尔可夫模型, fn (x0, x1, …, xn, xn+1) 为定义在Rn+2上有界Borel可测函数且│fn│≤M, 若{an, n≥0}是趋向于无穷的增序列, 且对任意自然数n,
│
则
证明 令
Zk=fk (X0, …, Xk, Yk) -E[fk (X0…, Xk, Yk) │Xk-1], k≥1 (2.3)
设Fn如引理1的定义, 下证{Zk, k≥1}是一个鞅差序列, 因为E[fk│Xk-1]为Fk-1可测的, 故
由引理1, (2.3) 式与 (2.4) 式, 有
E[Zk│Fk-1]=0 a.s. k≥1 (2.5)
故序列{Zk, k≥1}是一个鞅差序列.又
由条件期望Jensen的不等式, 有
故由 (2.1) 式与 (2.7) 式有
由 (2.1) 式有
由 (2.1) 式, (2.3) 式与 (2.8) 式, 有
由 (2.10) 式与引理2有
由 (2.3) 式与 (2.11) 式, 即得 (2.2) 式。
参考文献
[1]Leroux B G.Maximum-likelihood estimation for hidden Markov mod-els.Stoc hastic Processes and their Appl, 2002;40:127—143
[2]Bickel P J.Ritov Y, Ryden T.Aymptotic nomality of the maximum-like lihood estimator for general hidden Markov models.The Annals of Statistics, 2005;26 (4) :1614—1635
[3]Bates B C, Charles P, Hughes J P.Stochastic down-scaling of nu-merrical climate model simulations.Environmental Modeling soft-ware, 2006;13:325—331
[4]Lacruz B, Lasala P, Lekuona A.ynamic graphical models and nonhomo-geneous hidden Markov models.Stat Proba Letts, 2006;49:377—385
[5]Shaohua.Aconclusion about the transition matrix.World acdemic U-nion, 2006;112—115Theory, 1990;36 (5) :1006—1018
[6]Jin Shaohua.关于极限定理的一个结果及其推广.数学的实践与认识, 2007;37 (13) :118—123
[7]龚光鲁, 钱敏平.应用随机过程教程.北京:清华大学出版社, 2004:249
[8]Chung K L.A course in probability theory.New York;Academic Press, 1974
马尔可夫 篇8
随着信息化社会的不断发展, 人们对微软公司的Excel已不再陌生, 但是它的很多高级且强大的功能并不为一般用户所掌握。特别是与Excel绑定的VBA (Visual Basic For Application) , 它不仅可以实现程序自动化, 创建针对性强、实用性强和效率高的解决方案, 而且, 还可以将Excel用作开发平台, 实现更加复杂的应用程序系统开发[1]。在此利用语法简单、容易理解和掌握的VBA实现股票预测实例, 旨在普及和推广电脑的应用, 使程序设计和系统开发大众化, 让人们逐渐消除对 “编程”或“代码”的害怕心理, 能用VBA开发和解决本领域中的实际问题, 实现应用程序自动化, 提高办事效率, 减小工作量, 适应信息化社会不断发展的需要。
1马尔可夫预测法简介
1.1 马尔可夫过程和马尔可夫预测法
预测是人们根据历史资料和现实, 利用已经掌握的知识和手段, 对事物的未来或未知状况进行的事前推知或判断。其中马尔可夫法是预测中比较常用的方法。
马尔可夫 (Markov) 预测法, 简称马尔可夫法的基本原理是对事件的全面预测, 不仅要能够指出事件发生的各种可能结果, 而且还必须给出每一种结果出现的概率[2]。在事件的发展过程中, 若每次状态的转移都仅与前一时刻的状态有关, 而与过去的状态无关, 或者说状态转移过程是无后效性的, 则这样的状态转移过程就称为马尔可夫过程。所谓“无后效性”, 是指过去对未来无后效, 而不是指现在对未来无后效。马尔可夫链是与马尔可夫过程紧密相关的一个概念。马尔可夫链指出事物系统的状态由过去转变到现在, 再由现在转变到将来, 一环接一环, 像一根链条, 而作为马尔可夫链的动态系统将来是什么状态, 取什么值, 只与现在的状态、取值有关, 而与它以前的状态、取值无关[3,4]。因此, 运用马尔可夫链只需最近的或现在的动态资料便可预测将来。马尔可夫预测法就是应用马尔可夫链来预测未来变化状态。
马尔可夫预测法, 就是一种预测事件发生的概率的方法。它是基于马尔可夫链, 根据事件的目前状况预测其将来各个时刻 (或时期) 变动状况的一种预测方法。马尔可夫预测法是对地理、天气、市场等进行预测的基本方法。
1.2 转移概率和转移概率矩阵
转移概率就是用概率来描述事物状态转移的可能性大小, 常用“P ”表示[1]。
矩阵是指由许多个数组成的一个数表。每一个数称为矩阵的元素。其表示方法如下:
这是一个由m行n列的数构成的矩阵, Pij表示矩阵中的第i行与第j列交叉点上的元素, 矩阵中的行数与列数可以相等, 也可以不等。当它们相等时, 矩阵就是一个方阵。
由转移概率组成的矩阵就是转移概率矩阵。也就是说构成转移概率矩阵的元素是转移概率。转移概率矩阵有以下特征[5]: 矩阵中每个元素 (转移概率) 均为非负数;矩阵中每一行转移概率之和等于1。
1.3 几个基本概念
状态 某一事件在某个时刻 (或时期) 出现的某种结果。随着研究的事件及其预测的目标不同, 状态有不同的划分方式[6]。
状态转移过程 事件的发展, 从一种状态转变为另一种状态, 即为状态转移[6]。
一步转移概率 设系统有N个状态Ei (i=1, 2, …, N) , 以状态变量xt=i表示在时刻tn处于Ei (i=1, 2, …, N) , 如果系统在时刻tn处于Ei而在时刻tn+1转移到Ej的概率只与Ei有关而与tn以前处的状态无关, 则此概率可表示为:
Pij=P (Ei→Ej) =P ( xn+1=j∣xn=i) (0≤ Pij ≤1 , ∑ Pij=1) , 称Pij为一步转移概率。
预测模型[7]S (k+1) =S (k) P
起始状态概率分布 当预测模型中k=0时, 反映的是系统在初始时状态概率的分布情况, 称为起始状态概率分布。
k步转移概率 假设系统有N个状态Ei (i=1, 2, …, N) , 用Pi表示系统在k时期处于状态Ei (i=1, 2, …, N) 的概率, 所有概率所构成的向量, 称为状态概率向量。其中: 0≤Pi (k) ≤1 (i=1, 2, …, N) , ∑ Pi (k) =1, 由S (k+1) =S (k) P 可得递推关系S (k) =S (0) P×P×…×P (k个P相乘) , 称为k步转移概率[8]。
1.4 马尔柯夫预测法的步骤
根据上述分析, 马尔柯夫预测法的步骤应为[9,10,11]:
(1) 确定系统的状态Ei和S (0) ;
(2) 确定P;
(3) 进行预测:S (k) =S (0) Pk 。
假设系统的终极状态的状态概率为m1, m2, …, mN, 则: (m1, m2, …, mN) = (m1, m2, …, mN) P , 也就是说, 所谓终极状态就是再也不会改变的概率状态, 按照固定的概率进行转移。
2用VBA实现马尔柯夫预测法对股票预测
以股票收盘价状态 (上升、平盘、下降) 为对象进行预测, 选取前23天的行情作为样本来分析预测。其中前23天中上升9天, 平盘1天, 下降13天。因为第23天的状态是下降而无状态转移, 所以下降记为12天。其中由上升转移为上升的有2次, 由上升转移为平盘的有0次, 由上升转移为下降的有7次;平盘转移为上升的有0次, 由平盘转移为平盘的有0次, 由平盘转移为下降的有1次;由下降转移为上升的有6次, 由下降转移为平盘的有1次, 由下降转移为下降的有5次。由此, 可制作如图1所示。
假设想知道3天后股票的情况, 要先得出一步概率转移矩阵, 其程序代码如下:
Sub OneStep ()
其结果如图2所示。
因为要看3天后股票的情况, 所以需要做一个通式, 只要输入天数, 不论是几天后的情况都能得到。该例先用下列代码把天数限制在一周内:
然后, 用下列代码来计算矩阵的k次幂, 即k步转移概率:
该代码的作用是只要每次C15单元格数字发生变化, 工作表便会调用函数AnyStep () , 从而得出k步转换概率矩阵, 得到相关月份的所有数据。如在C15单元格中输入3, 按回车或者将鼠标放置其它处, 便能得到第3天的预测情况, 如图3所示。
本例的天数被限制在一周内, 如果在C15单元格输入非1～7的数字, 则会报告错误。提示输入合法数值, 如图4所示。
3结语
通过应用证明, 该程序比较适用、操作简单、预测结果和实际结果有很好的一致性。当然, Excel在预测中的应用是非常广泛的, 如组织内部人力资源供给的预测、市场占有率的预测和销售期望利润的预测等。基于预测的经常性, 方法的固定性, 能用Excel的时候不妨都用VBA, 尽管第一次可能稍微多花些时间, 但却会极大地减少以后的工作量。另外, 由于VBA编辑器中提供了丰富的控件和完备的语言系统, 能实现较大的管理、控制系统, 用户可以根据需要创建自己的VBA应用程序, 所以对于非专业软件开发人员, 用Excel及VBA开发本领域的系统是一个非常不错的选择。
参考文献
[1]韩泽坤, 朱瑞亮.Excel VBA高效办公范例应用[M].北京:中国青年出版社, 2004.
[2]张宗震.马尔可夫预测法基本原理解析与应用[J].成都电子机械高等专科学校学报, 2001 (2) :21-26, 31.
[3]倪杰.吸收态马尔可夫链及其应用[J].商业研究, 2000 (7) :75-76.
[4]黄勇, 吴庆文.建筑陶瓷出口营销预测的马尔柯夫链模型[J].中国陶瓷, 2004 (4) :8-11.
[5]郝海, 张双德.马尔可夫模型及应用[J].中国民航学院学报, 2003 (Z1) :11-13.
[6]Trinkl F A.Stochastic Analysis of Programs for MentallyRetarded[M].Operations Research, 1984.
[7]Meredith J.A Markovian Analysis of Geriatric Ward[M].Management Science, 1979.
[8]Barro R J, Salad-Martin X.Economic Growth.1995.
[9]张少文.黄河流域天然年径流变化特性分析及其预测[D].成都:四川大学, 2005.
[10]姜冠杰, 黄敏, 金天昊.基于马尔可夫过程的虚拟企业动态风险评价方法的研究[J].鞍山科技大学学报, 2005 (Z1) :12-15.