Markov链预测法资料.pdf

《Markov链预测法资料.pdf》由会员分享，可在线阅读，更多相关《Markov链预测法资料.pdf（21页珍藏版）》请在三一文库上搜索。

1、1 我们仔细阅读了中国大学生数学建模竞赛的竞赛规则. 我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮 件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问 题。 我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他 公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正 文引用处和参考文献中明确列出。 我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反 竞赛规则的行为，我们将受到严肃处理。 我们参赛选择的题号是（从A/B/C/D 中选择一项填写） ： B 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学

2、校（请填写完整的全名）： 贵州民族学院 参赛队员(打印并签名 ) ：1. 龚道杰 2. 张凤 3. 姚肖伟 指导教师或指导教师组负责人(打印并签名 )： 日期： 2009 年 7 月 25 日 年凝冻日数的 Markov 链预测法 4# 【摘要】 本文根据所给数据，利用 Markov链建立了预测年凝冻日数的模型，分别从 整体和局部两个角度进行分析。 2 首先，我们直接以年凝冻日数为依据，对其进行K-均值聚类分析，划分 状态。用频率估计概率的方法， 估算出一步转移概率矩阵， 1/ 65/ 6 5/ 3328/ 33 P， 然后建立 Markov 链模型 ( ) 1/ 65/6 ( )(0)(0)

3、 5/ 3328 /33 n n P nPPP。以 2008 年 作 为 初 始 状 态 ， 估 计 出2009年 凝 冻 日 数 所 处 状 态 为 (1)(0)PPP0.1520.848 。按 K-均值标准可知， 即 2009 年凝冻的天数在 15 天以内的可能性为84.8%，在 15 天以上的可能性为15.2%。 由于上述模型选取的是以年为单位的数据，只能估计出2009年的凝冻日 数所处区间。为提高精度，我们选取2000-2008年的具体凝冻天数和日期，记每 一天只存在两种状态，出现雨凇为状态1，否则为状态 0。然后由相邻两年间的 状态转移变化，得出一步转移概率矩阵 i P ，1,2,.

4、,8i。由这 8个一步转移概率 矩阵，根据一步转移矩阵P的 n次方与 n 步转移概率矩阵 ( )n P之差的范数和达到 最小的准则，选出优化后的一步转移概率矩阵 0.95000.0500 * 0.78890.2111 P , 再次建立 Markov 链模型。以 2008 年为初始状态，预测 2009年的概率分布为 * (2009)(2008)0.91060.0894PPP，由频率稳定于概率，知2009年凝冻 天数的估计值为14天。 关键词：Markov 链转移概率矩阵频率估计概率 1. 问题提出 1.1 背景知识 凝冻是指冬季出现的温度低于0有过冷却降水或固体降水和结冰现象 发生的天气现象，即

5、气象台所说的出现雨凇的天气。雨凇的形成与气温，降水 量，湿度等因素有关，超冷却的降水碰到温度等于或低于零摄氏度的物体表面 使所形成玻璃状的透明或无光泽的表面粗糙并覆盖层，就叫做雨凇。其造成的 危害巨大，高压线塔的倒塌，电力瘫痪，交通瘫痪，农作物的冻亡等。因 而对出现雨凇天气的预测显得尤为重要。 3 1.2 问题分析 根据所给 1969-2008 年的数据，建立一个年凝冻日数的预测模型， 预测 2009 年的凝冻日数，并作出误差分析。数据给出了是否出现雨凇与气 温、降水量、湿度、气压和风速的关系，而雨凇的出现是一个随机过程， 与多个因素有关，且受干预变量的影响，因而传统的回归分析方法，效果 不好

6、，而 Markov 链构造模型不需要从复杂的预测因子中寻找各因素之间的 相互规律，只需要考虑事件本身的演变特点，通过计算转移概率矩阵来预 测内部状态的变化。 2. 建模准备 2.1 数据分析与处理 以年为单位，统计出现雨凇的天数，见表1： 年份日数年份日数年份日数 1969 15 1983 3 1997 5 1970 6 1984 27 1998 4 1971 0 1985 0 1999 0 1972 8 1986 3 2000 10 1973 1 1987 3 2001 6 1974 20 1988 12 2002 9 1975 8 1989 8 2003 8 1976 8 1990 6 2

7、004 17 1977 16 1991 0 2005 18 1978 6 1992 1 2006 10 1979 4 1993 6 2007 8 1980 8 1994 3 2008 37 1981 10 1995 0 1982 7 1996 8 2.2 Markov 链预测的理论基础 2.2.1 Markov链定义 （Markov 链） 1 随机过程 Xn，0,1,2.n称为 Markov 链，若它只取有 限或可列个值 012 ,.EE E（我们以0,1,2,.来标记 01 ,.,E E并称它们是过程 的状态，0,1,2.或者其子集记为S，称为过程的状态空间） . 对任意的0n 及状态 01

8、1 , , , ., n i j iii有 4 1 n P Xj 00112211 ,., nnn XiXiXiXiXi = 1 n P Xj . n Xi（5.1.1 ） 式（5.1.1 ）刻画了 Markov 链的特性，称为 Markov 性。 2.2.2 转移概率矩阵 由转移概率组成的矩阵，形如 000102 101112 202122 . . () . ij ppp ppp Pp ppp 称 P为转移概率矩阵。且 ij p( ,)i jS有性质： (1)0, ,; ij pi jS (2)1,. ij j s piS【2】 2.2.3 （C-K方程） 对一切,0, ,n mi jS有

9、() (1); m nmn ijikkj ks pp p (1)(2) (2). nnn PP PP P PP (n) 其证明如下： () 0 | m n ijm n pp Xj Xi = 0 0 , m n p Xj Xi p Xi = 0 0 , m nm ks p Xj Xk Xi p Xi （全概率公式） = 00 00 , , m nmm ks m p Xj Xk Xip Xk Xi p Xip Xk Xi = 00 |,| m nmm ks p Xj Xk Xi p Xk Xi = ( )()nm kjik ks pp = ()( )mn ikkj ks pp【3】 5 2.3.4

10、 传统的频率估计概率估算一步转移概率矩阵的方法为： 已知系统存在 n 种状态，状态空间为S=0,1,2,n. 假设在次观 测中，系统处于第i种状态共有 i n 次，显然 1 n j j Nn 用 ij n 表示系统从状态i 经过一步转移到状态j的频数，显然有 1 ,( ,), n iijij j nni jS n 组成的矩阵 () ij n 称为转移频数矩阵。 将转移频数矩阵的第i行第j列元素 ij n 除以i行各元素总 和所得的值称为转移概率，记为, , ij p i jI 。即有/ ijiji pnn ，于是我们得到 用频率估计出一步转移概率矩阵P【】 3. 符号说明 符号说明 i t第j

11、期的概率分布 ij P从状态i到状态j的转移概率 I状态空间且0,1I i t频率 4 模型的建立 4.1 模型假设 1）雨凇的年出现次数是一簇依赖于时间的随机变量，其变化过程是一个随 机过程； 2 ）该随机过程具有无后效性； 3 ）雨凇年出现次数状态的一步转移概率矩阵只与时间差有关，与时间起点 无关。 4.2 模型建立 4.2.1 以表 1 为基础，建立 Markov 链预测模型： 6 1 ）利用 SPSS 软件，以 K均值聚类法将过去的年凝冻日数分为2 个区间， 确定每年凝冻日数的状态, 见表 2： 2）根据表 2，以频率估计概率的方法，计算一步转移概率矩阵。 出现状态1的次数为716，出

12、现状态2的次数为33。由1转为1的 次数为1，故转移概率 11 1/ 6P；由1转为2的次数为5，故转移概率 12 5/ 6P； 由2转为1的次数为5，故转移概率 21 5/ 33P；由2转为2的次数为28，故转 移概率 22=28/33 P。 由此可得雨凇年出现次数状态的一步转移概率矩阵为： 1/ 65/ 6 5/ 3328/ 33 P； Markov 链的基本原理就是利用初始状态概率向量和状态转移概率矩 阵来推知预测对象将来一个时期所处的状态。 记 0 (0), j PP XjjS，则有 12 (0)(0),(0),.,(0),.) j PPPP，称它为 Markov 链的初始分布，显然有

13、(0)1 j jS P。由上述 C-K 方程 可知 Markov 链 在任一时刻n的一维分布由初始分布(0)P和n步转移概率矩阵所确定。 即 Markov 链的预测模型为 ( ) 1/ 65/6 ( )(0)(0) 5/3328/33 n n P nPPP。 （1） 4.2.2 根据所建 Markov 链模型，进行预测 用 2008 年凝冻天数作为初始状态，即(0)01P. 利用模型（ 1） 式，计算可得 2009年凝冻天数的一维分布为： 7 1/ 65/ 6 (1)(0)01 5/ 3328/ 33 PPP0.1520.848 这表明 2009 年的凝冻天数所处的状态为1 的概率为 0.15

14、2，状态为 2 的概率为 0.848. 由之前 SPSS 软件的 K-均值聚类可知，凝冻的天数在15 天以 内的可能性为 84.8%，在 15天以上的可能性为15.2%。 4.3 模型检验和结果分析 该模型虽然预测出了2009 年凝冻日数的范围， 并计算出其以 84.8%的概 率稳定于该状态，却无法的估计出2009 年凝冻的具体天数。 由于凝冻基本发生在1 月、2 月、3 月、11月、12 月，而 2009 年前三 个月的历史天气数据可以查得，见数据1 【5】 2009 年贵阳雨淞出现的次数 日期雨淞出现日期雨淞出现日期雨淞出现 1-1 1 2-1 0 3-1 0 1-2 0 2-2 0 3-

15、2 0 1-3 0 2-3 0 3-3 0 1-4 1 2-4 0 3-4 0 1-5 0 2-5 0 3-5 0 1-6 1 2-6 0 3-6 0 1-7 1 2-7 0 3-7 0 1-8 0 2-8 0 3-8 0 1-9 0 2-9 0 3-9 0 1-10 0 2-10 0 3-10 0 1-11 0 2-11 0 3-11 0 1-12 0 2-12 0 3-12 0 1-13 0 2-13 0 3-13 0 1-14 0 2-14 0 3-14 0 1-15 0 2-15 0 3-15 0 1-16 0 2-16 0 3-16 0 1-17 0 2-17 0 3-17 0 1-

16、18 0 2-18 0 3-18 0 1-19 0 2-19 0 3-19 0 1-20 0 2-20 0 3-20 0 1-21 0 2-21 0 3-21 0 1-22 0 2-22 0 3-22 0 1-23 0 2-23 0 3-23 0 1-24 1 2-24 0 3-24 0 1-25 1 2-25 0 3-25 0 1-26 1 2-26 0 3-26 0 1-27 1 2-27 0 3-27 0 8 1-28 0 2-28 0 3-28 0 1-29 0 3-29 0 1-30 0 3-30 0 1-31 0 3-31 0 由数据 1 可得， 2009 年发生凝冻 1 月天数为

17、 8 天，2 月天数为 0 天，3 月天数为 0 天。 对题目所附数据做简单统计分析，见表 3 表 3 40以来各月份出现凝冻的天数 1 月份2 月份3 月份11 月份12 月份 雨凇出现日 数总和 16410811343 雨凇出现日 数平均 数 4.12.70.2750.0751.075 根 据 上 表 可 知 ， 凝 冻 发 生 在 前 三 个 月 的 频 率 为 1 (164 108 11)/3290.8602t ,发生 在 后 两 个 月的 频 率 为 2 (343)/3290.1398t。即凝冻发生在11 月、12 月的天数和远小于1 月、 2 月、 3 月的天数和，粗略估计 200

18、9年 11 月、 12 月的天数和小于 8 天， 则 2009 全年凝冻天数小于15 天。与模型（ 1）非常吻合。 Markov 链预测模型成功的关键在于转移概率矩阵的可靠性，因此模型的构 造需要足够多的准确的统计数据，而本题提供了40 个年度凝冻日数的数据偏少， 会影响预测精度。本题在求转移概率矩阵的时候，采用的是传统的估算方法，先 假设已知随机过程在n 种状态的观测次数及系统从当前时刻向下一时刻转移次数 的情况下，用频率估计概率的方法估算出一步转移概率矩阵。但在实际情况下， 没有足够的观测次数，会导致一步转移概率矩阵和真实值相差很大。 对于本题，如果改为从具体每天是否出现雨凇的状态考虑，4

19、0 年的海量数 据，将会极大提高我们模型的估计精度。 5模型的改进 5.1 数据分析 用 SAS软件对表 1 的数据作时序图和自相关图（程序见附件），检验其 平稳性。 9 时序图： 自相关图： 结合时序图与自相关图分析，以年为单位， 其凝冻天数基本上平稳。我 们不妨取 20002008年的数据进行分析预测。 5.2 对数据处理并再次建模 从 2000-2008 年，2 月有 28 或 29 天，为计算方便，统一取28 天，则每 年关于凝冻的数据有 151 个。 10 2000 年出现雨凇的天数 1 月2 月3 月11月12 月 00000 01000 01000 01000 01000 000

20、00 00000 00000 00000 00000 00000 00000 00000 00000 00000 10000 10000 10000 00000 00000 00000 00000 00000 00000 00000 10000 00000 10000 1000 0000 000 2001 年出现雨凇的天数 1 月2 月3 月11月12月 00000 00000 00000 11 00000 00000 00000 00000 00000 10000 00000 00000 00000 00000 00001 10000 00000 00000 00001 00000 0000

21、0 00000 00000 00000 00000 10000 10000 00000 00000 0000 0000 000 观察上述两组数据， 对于每一天雨凇的出现只存在两种情况，出现或否 , 即两种状态0,1I. 现在我们以天为单位，以2000-2001 年的数据为例进行 分析。2000 年的统计数据有151 个，其中处于状态0的有 141 个，处于状态 1 的有 10 个。 记其分布为(2000)(141/151 10/151)P； 同理可知(2001)(145/1516/151)P； 从 2000 年转移至 2001 年， 状态00有 138个；状态01有 3 个； 状态10有 9个

22、；状态11有 1 个。 12 则 2000 年到 2001年的转移概率矩阵为 1 138/1413/141 9/101/10 P； 既有 1 138/1413/141 (2000)*141/151 10/151(2001) 9/101/10 PPP。 同理 2001 年到 2002年转移概率矩阵为 2 138/1468/146 4/51/ 5 p； 2002年到 2003 年转移概率矩阵为 3 134/1428/142 10 p； 2003年到 2004 年转移概率矩阵为 4 127/14316/143 7/81/8 p； 2004年到 2005 年转移概率矩阵为 5 119/13415/13

23、4 14/173/17 p； 2005年到 2006 年转移概率矩阵为 6 127/1336/133 7/ 92/ 9 p； 2006年到 2007 年转移概率矩阵为 7 133/1418/141 10 p； 2007年到 2008 年转移概率矩阵为 8 111/14332/143 3/85/8 p。 由以上八个转移概率矩阵可以看出，实际生活中相邻时刻的一步转移概率 矩阵并不是完全相等的，为了能得出一个尽量精确的一步转移概率矩阵来预测 2009年的数据，我们需要对上述转移概率矩阵进行优化。余波等人给出了利用 最优化的思想，使一步转移矩阵P的 n次方与 n步转移概率矩阵 ( )n P之差的范数

24、13 和达到最小的准则，建立模型如下： 目标函数： 1 min()| n i i f PPP（2） 约束条件： 1 1,1,2,., 0, ,1,2,., n ij j ij pin Pi jn 。 【6】 （矩阵范数）对任意, m n A BR称| |为 m n R空间的矩阵范数，指 |满足： (1)| 0;| 00 (2) | | (3)| | AAA AA ABAB 对任意C 设() ij AaM ， 定义 1/2 2 ,1 | n ij i j aA。【7】 结合本题数据，利用模型（2）建立规划求解： 优化后的一步转移概率矩阵记为 * 0. 9 5 0 00. 0 5 0 0 ( 2

25、0 0 9 )( 2 0 0 8 ) 1 1 4 / 1 5 13 7 / 1 5 10. 9 1 0 60. 0 8 9 0. 7 8 8 90. 2 1 1 1 PPP 目标函数： 8 * 1 m i n()| | | i i fPPP 约束条件： 由 MATLAB 软件求解得 0. 9 5 0 00. 0 5 0 * 0. 7 8 8 90. 2 1 1 1 P ; 以 2008 年基期，预测 2009 年的概率分布有 * ( 2 0 0 9 )( 2 0 0 8 )0. 9 1 0 60. 0 8 9 4PPP1 5 1 *( 2 0 0 9 )1 3 7. 4 9 4 21 3. 5

26、 0 5 8P 则预测 2009 年发生凝冻的天数为 14 天。 6、模型的优缺点分析 本文从整体和局部两个角度分别建模，两个模型的结果大致吻合，说明 Markov 链模拟的效果还不错。 对于本题，雨凇的出现可以看成是一个随机过程， 而影响雨凇出现的因素太多，若直接分析影响因素，会十分麻烦，且由于干预 变量会使模型的精度大为降低。 Markov 链模型的优点在于它不需要从复杂的预测因子中寻找个因素之间 的相互规律，直接通过概率矩阵来预测内部状态的变化。可是Markov 链对转移 概率矩阵的要求很高，实际中，由于观测数据的限制，误差影响，会造成转移 14 概率矩阵精度的降低。因此为保证Marko

27、v 链模拟的准确性，需要收集足够多的 准确的数据。 参考文献 【1】 张波，张景肖 . 应用随机过程 . 北京：清华大学出版社，2008 年，P74 【2】 张波，张景肖 . 应用随机过程 . 北京：清华大学出版社，2008 年，P75. 【3】 张波，张景肖 . 应用随机过程 . 北京：清华大学出版社，2008 年，P81. 【4】 于波, 陈希镇 , 华栋. 马尔柯夫链在农作物年景预测中的应用. 统计与决 策 2007 年第 21期. 【5】 历史天气数据http:/ 2009 年 7 月 23日. 【6】 于波, 陈希镇 , 华栋. 马尔柯夫链在农作物年景预测中的应用. 统计与决 策 2

28、007 年第 21期. 【7】 ppt Chapter1.2 向 量 范 数 与 矩 阵 范 数 武 汉 大 学 精 品 课 程 http:/ 年 7 月 23 日 附录 1、 程序 （时序图） data text1; input days; time=intnx(year, 1969,_n_-1); cards ; 15 6 0 8 1 20 8 8 16 6 4 8 10 7 3 27 0 3 3 12 8 6 0 1 6 3 0 8 5 4 0 10 6 9 8 17 18 108 37 ; procgplotdata =text1; plot days*time; symbolc=bl

29、ack v=star i =join; run ; （自相关图） data text2; input freq; time=intnx(year, 1969,_n_-1); cards ; 15 6 0 8 1 20 8 8 16 6 4 8 10 7 3 27 0 3 3 12 8 6 0 1 6 3 0 8 5 4 0 10 6 9 8 17 18 10 8 37 15 ; procarimadata =text2; identifyvar =freq; run ; 2. 数据 2002年出现雨凇的天数 1 月2 月3 月11月12 月 00000 00000 00000 00000 00

30、000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 10000 10000 00000 00000 00000 10000 10001 00001 00001 00001 0000 0001 000 16 2003年出现雨凇的天数 1 月2 月3 月11月12 月 00000 01000 00000 00000 10100 10100 00000 00000 10000 00000 00000 01000 00000 00000 00000 00000 00000 00000 000

31、00 00000 00000 00000 00000 00000 00000 00000 10000 00000 0000 0000 000 2004年出现雨凇的天数 1 月2 月3 月11月12 月 00000 01000 01000 17 01000 01000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 10000 00000 10000 00000 10000 10000 00001 10001 10000 00001 00001 0001 0001 001 2005 年出现雨凇的天数

32、1 月2 月3 月11月12 月 11000 00000 00000 00000 00001 00001 00000 00000 10000 11000 18 11000 01000 01000 00000 00000 00000 01000 00000 00000 10000 10000 10000 00000 00000 00000 10000 10000 00000 0000 0000 100 2006 年出现雨凇的天数 1 月2 月3 月11月12 月 00000 00000 00000 01000 10000 10000 10000 00000 00000 00000 00000 0

33、0000 00100 00000 00000 00000 01000 19 01000 00000 10000 10000 10000 00000 00000 00000 00000 00000 00000 0000 0000 000 2007 年出现雨凇的天数 1 月2 月3 月11月12 月 00000 00000 10000 00000 00000 00000 00100 00000 00000 00000 00000 10000 10000 10000 10000 10000 10000 00000 00000 00000 00000 00000 00000 20 00000 00000 00000 00000 00000 0000 000 000 2008 年出现雨凇的天数 1 月2 月3 月11月12 月 01000 01000 01000 01000 01000 01000 01000 01000 01000 01000 01000 01000 11000 10000 10000 10000 10000 10000 10000 10000 10000 10001 10001 10001 10000 11000 10000 10000 1000 1001 21 100