我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的公开的资料（包括网上查到的资料）文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。
, 如果引用别人的成果或其他
，必须按照规定的参考文献的表述方式在正
我们参赛选择的题号是（从A/B/C/D中选择一项填写）：B
：
我们的参赛报名号为（如果赛区设置报名号的话）所属学校（请填写完整的全名）参赛队员
(打印并签名) ：1.
2.
：贵州民族学院龚道杰张凤
3. 姚肖伟
指导教师或指导教师组负责人
(打印并签名)：
日期： 2009
年 7
月 25
日
年凝冻日数的Markov链预测法4#
【摘要】
本文根据所给数据，利用Markov链建立了预测年凝冻日数的模型，分别从整体和局部两个角度进行分析。
1
首先，我们直接以年凝冻日数为依据，对其进行K-均值聚类分析，划分
1/6
5/6
5/3328/33
n
状态。用频率估计概率的方法，估算出一步转移概率矩阵，P
1/65/33
5/628/33
，
然后建立Markov链模型P(n)
P(0)P
(n)
P(0)
。以2008年
作为初始状态，估计出
P(1)
P(0)P
2009年凝冻日数所处状态为
0.1520.848。按K-均值标准可知，即2009年凝冻的天数在
15天以内的可能性为84.8%，在15天以上的可能性为15.2%。
由于上述模型选取的是以年为单位的数据，只能估计出2009年的凝冻日
数所处区间。为提高精度，我们选取2000-2008年的具体凝冻天数和日期，记每一天只存在两种状态，出现雨凇为状态1，否则为状态0。然后由相邻两年间的状态转移变化，得出一步转移概率矩阵矩阵，根据一步转移矩阵
Pi，i
1,2,...,8。由这8个一步转移概率
(n)
P的n次方与n步转移概率矩阵P
之差的范数和达到
最小的准则，选出优化后的一步转移概率矩阵P*
0.95000.05000.78890.2111
,再次建立Markov链模型。以2008年为初始状态，预测
2009年的概率分布为P(2009)
P(2008)P
*
0.91060.0894，由频率稳定于概率，知2009年凝冻
天数的估计值为14天。
关键词：Markov链转移概率矩阵频率估计概率
1. 问题提出
1.1背景知识
凝冻是指冬季出现的温度低于0℃有过冷却降水或固体降水和结冰现象
发生的天气现象，即气象台所说的出现雨凇的天气。雨凇的形成与气温，降水量，湿度等因素有关，超冷却的降水碰到温度等于或低于零摄氏度的物体表面使所形成玻璃状的透明或无光泽的表面粗糙并覆盖层，就叫做雨凇。其造成的危害巨大，高压线塔的倒塌，电力瘫痪，交通瘫痪，农作物的冻亡等。因而对出现雨凇天气的预测显得尤为重要。
2
1.2问题分析
根据所给1969-2008年的数据，建立一个年凝冻日数的预测模型，
预测2009年的凝冻日数，并作出误差分析。数据给出了是否出现雨凇与气温、降水量、湿度、气压和风速的关系，而雨凇的出现是一个随机过程，与多个因素有关，且受干预变量的影响，因而传统的回归分析方法，效果不好，而Markov链构造模型不需要从复杂的预测因子中寻找各因素之间的相互规律，只需要考虑事件本身的演变特点，通过计算转移概率矩阵来预测内部状态的变化。
2. 建模准备
2.1数据分析与处理
以年为单位，统计出现雨凇的天数，见表年份1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982
日数15 6 0 8 1 20 8 8 16 6 4 8 10 7
年份1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996
日数3 27 0 3 3 12 8 6 0 1 6 3 0 8 1：
年份1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008
日数5 4 0 10 6 9 8 17 18 10 8 37
2.2 Markov链预测的理论基础2.2.1 Markov链定义
（Markov链）
[1]
随机过程{Xn，n0,1,2...｝称为Markov链，若它只取有
限或可列个值E0,E1,E2,...（我们以｛0,1,2,...｝来标记E0,E1,...,并称它们是过程的状态，｛0,1,2...｝或者其子集记为S，称为过程的状态空间）.对任意的n及状态i,j,i0,i1...,in1,有
0
3
P{Xn
=P{Xn
1
j︱X0j︱Xn
i0,X1i}.
i1,X2i2,...,Xn
1
in1,Xn
i｝
1
（5.1.1）
式（5.1.1）刻画了Markov链的特性，称为Markov性。2.2.2 转移概率矩阵
由转移概率组成的矩阵，形如
p00
P
(pij)
p10p20
p01p11p21
p02p12p22
.........
称P为转移概率矩阵。且pij(i,j
(1)pij
S)有性质：
0,i,jpij
1,i
S;S.
【2】
S有
(2)
js
2.2.3（C-K方程）对一切n,m0,i,j
(1)p
(mn)ij
ks
pp;PP
(n1)
miknkj
(2)P
(n)
PPP
(n2)
P.
n
其证明如下：
p
=
(mn)ij
pXmj,X0
ij,XmpX0
i
n
j|X0
i
pXm
n
pX0pXm
ks
n
=
k,X0ik,X0ik,X0
i
（全概率公式）ipXm
pXmipXm
k,X0k,X0k|X0
iii
=
ks
pXm
n
j,XmpX0
=
ks
pXmp
ks
(n)kj
n
j|Xm
= =
ks
pp
(m)ik
p
(m)ik(n)kj
【3】
4
2.3.4 传统的频率估计概率估算一步转移概率矩阵的方法为：
已知系统存在n种状态，状态空间为S={0,1,2,…n}.假设在Ｎ次观
n
测中，系统处于第i种状态共有ni次，显然N
j1n
nj．用nij表示系统从状态i
经过一步转移到状态
j的频数，显然有ni
j1
nij,(i,j
S),nij组成的矩阵(nij)
称为转移频数矩阵。将转移频数矩阵的第i行第j列元素nij除以i行各元素总和所得的值称为转移概率，记为用频率估计出一步转移概率矩阵
pij,i,j
I。即有pij
nij/ni，于是我们得到
P．【４】
3.符号说明
符号
说明
第j期的概率分布
Pij
I
ti
从状态i到状态j的转移概率状态空间且I频率
{0,1}
ti
4．模型的建立
4.1 模型假设
1）雨凇的年出现次数是一簇依赖于时间的随机变量，其变化过程是一个随机过程；
2）该随机过程具有无后效性；
3）雨凇年出现次数状态的一步转移概率矩阵只与时间差有关，与时间起点无关。4.2模型建立
4.2.1以表1为基础，建立Markov链预测模型：
5
1）利用SPSS软件，以K均值聚类法将过去的年凝冻日数分为确定每年凝冻日数的状态,见表2：
2个区间，
2）根据表2，以频率估计概率的方法，计算一步转移概率矩阵。
出现状态1的次数为71
次数为1，故转移概率P11
6，出现状态2的次数为33。由1转为1的
1/6；由1转为2的次数为5，故转移概率P125/6；
由2转为1的次数为5，故转移概率P21移概率P22=28/33。
5/33；由2转为2的次数为28，故转
由此可得雨凇年出现次数状态的一步转移概率矩阵为：
P
1/6
5/6
5/3328/33
；
Markov链的基本原理就是利用初始状态概率向量和状态转移概率矩
阵来推知预测对象将来一个时期所处的状态。
记Pj(0)
P{X0
j},j
S，则有P(0)
(P1(0),P2(0),...,Pj(0),...)，称它为
Markov链的初始分布，显然有
jS
Pj(0)1。由上述 C-K方程可知Markov链
P(0)和n步转移概率矩阵所确定。
在任一时刻n的一维分布由初始分布
即Markov链的预测模型为
（1）
P(n)P(0)P
(n)
P(0)
1/65/33
5/628/33
n
。
4.2.2 根据所建Markov链模型，进行预测
用2008年凝冻天数作为初始状态，即
式，计算可得2009年凝冻天数的一维分布为：
6
P(0)01.利用模型（1）
P(1)P(0)P01
1/65/6
5/3328/33
0.1520.848
这表明2009年的凝冻天数所处的状态为1的概率为0.152，状态为2
15天以
的概率为0.848.由之前SPSS软件的K-均值聚类可知，凝冻的天数在内的可能性为84.8%，在15天以上的可能性为15.2%。4.3 模型检验和结果分析
该模型虽然预测出了2009年凝冻日数的范围，并计算出其以84.8%的概率稳定于该状态，却无法的估计出
2009年凝冻的具体天数。
1 【5】
由于凝冻基本发生在1月、2月、3月、11月、12月，而2009年前三个月的历史天气数据可以查得，见数据
日期
1-1 1-2 1-3 1-4 1-5 1-6 1-7 1-8 1-9 1-10 1-11 1-12 1-13 1-14 1-15 1-16 1-17 1-18 1-19 1-20 1-21 1-22 1-23 1-24 1-25 1-26 1-27
2009年贵阳雨淞出现的次数雨淞出现日期雨淞出现日期雨淞出现1 2-1 0 3-1 0 0 2-2 0 3-2 0 0 2-3 0 3-3 0 1 2-4 0 3-4 0 0 2-5 0 3-5 0 1 2-6 0 3-6 0 1 2-7 0 3-7 0 0 2-8 0 3-8 0 0 2-9 0 3-9 0 0 2-10 0 3-10 0 0 2-11 0 3-11 0 0 2-12 0 3-12 0 0 2-13 0 3-13 0 0 2-14 0 3-14 0 0 2-15 0 3-15 0 0 2-16 0 3-16 0 0 2-17 0 3-17 0 0 2-18 0 3-18 0 0 2-19 0 3-19 0 0 2-20 0 3-20 0 0 2-21 0 3-21 0 0 2-22 0 3-22 0 0 2-23 0 3-23 0 1 2-24 0 3-24 0 1 2-25 0 3-25 0 1 2-26 0 3-26 0 1 2-27 0 3-27 0
7
1-28 1-29 1-30 1-31 0 0 0 0
2-28 0
3-28 3-29 3-30 3-31 0 0 0 0
由数据1可得，2009年发生凝冻1月天数为8天，2月天数为0天，3月天数为0天。
对题目所附数据做简单统计分析，见
表3 401月份
雨凇出现日
数总和雨凇出现日
数平均数
根据上表可知，凝冻发生在前三个月的频率为t1t2
(16410811)/329(343)/329
0.8602 ,
发生在后两个月的频率为
4.1
2.7
0.275
0.075
1.075
164
2月份108
表3 3月份11
11月份
3
12月份
43
以来各月份出现凝冻的天数
0.1398。即凝冻发生在11月、12月的天数和远小于1月、
2月、3月的天数和，粗略估计2009年11月、12月的天数和小于8天，则2009全年凝冻天数小于15天。与模型（1）非常吻合。
Markov链预测模型成功的关键在于转移概率矩阵的可靠性，因此模型的构造需要足够多的准确的统计数据，而本题提供了假设已知随机过程在
40个年度凝冻日数的数据偏少，
会影响预测精度。本题在求转移概率矩阵的时候，采用的是传统的估算方法，先
n种状态的观测次数及系统从当前时刻向下一时刻转移次数
的情况下，用频率估计概率的方法估算出一步转移概率矩阵。但在实际情况下，没有足够的观测次数，会导致一步转移概率矩阵和真实值相差很大。
对于本题，如果改为从具体每天是否出现雨凇的状态考虑，据，将会极大提高我们模型的估计精度。
5．模型的改进
5.1数据分析
用SAS软件对表1的数据作时序图和自相关图（程序见附件），检验其
平稳性。
40年的海量数
8
时序图：
自相关图：
结合时序图与自相关图分析，以年为单位，其凝冻天数基本上平稳。我
们不妨取2000—2008年的数据进行分析预测。5.2 对数据处理并再次建模
从2000-2008年，2月有28或29天，为计算方便，统一取年关于凝冻的数据有 151个。
9
28天，则每
1月2000年出现雨凇的天数2月3月11月12月
00001001001001000000000000000000000000000000010010010000000000000000000000010000010010000
0
2001年出现雨凇的天数
1月
2月
3月
11月
0000000000000000000000000000000000000
0012月
000
10
0000000000000000000000000000000
000
0000010000010000000001100000
0000000000000000000000000
0000000000000000000000000000
000000000000000000000000000
0000000000100010000000000000
出现或否,
观察上述两组数据，对于每一天雨凇的出现只存在两种情况，即两种状态I
{0,1}.现在我们以天为单位，以
2000-2001年的数据为例进行
分析。2000年的统计数据有151个，其中处于状态0的有141个，处于状态1的有10个。
记其分布为同理可知
P(2000)
(141/15110/151)
；；
P(2001)(145/1516/151)
从2000年转移至2001年，
状态0状态1
00
有138个；有 9个；
状态0状态1
11
有3个；有1个。
11
则2000年到2001年的转移概率矩阵为
P1
138/1413/1419/10
1/10
；
138/1413/1419/10
1/10
既有P(2000)*P1
141/15110/151P(2001)。
同理 2001年到2002年转移概率矩阵为
p2
2002
138/1464/5
8/1461/5
；
年到2003年转移概率矩阵为
p3
134/1428/142
1
0
；
2003年到2004年转移概率矩阵为
p4
127/14316/1437/8
1/8
；
2004年到2005年转移概率矩阵为
p5
119/13415/13414/17
3/17
；
2005年到2006年转移概率矩阵为
p6
127/1336/1337/9
2/9
；
2006年到2007年转移概率矩阵为
p7
133/1418/141
1
0
；
2007年到2008年转移概率矩阵为
p8
111/14332/1433/8
5/8
。
由以上八个转移概率矩阵可以看出，实际生活中相邻时刻的一步转移概率
矩阵并不是完全相等的，为了能得出一个尽量精确的一步转移概率矩阵来预测2009年的数据，我们需要对上述转移概率矩阵进行优化。余波等人给出了利用最优化的思想，使一步转移矩阵
P的n次方与n步转移概率矩阵P
(n)
之差的范数
12
和达到最小的准则，建立模型如下：
n
目标函数：
minf(P)
||Pi
P||
（2）
i1
n
约束条件：
pij
1,i1,2,...,nj1
。【6】
Pij
0,i,j
1,2,...,n
（矩阵范数）对任意A,B
R
mn
,称||||为R
mn
空间的矩阵范数，指||||满足：
(1)||A||0;||A||0
A0
(2)||A||||||A||对任意
C
(3)||AB||||A||||B||
n
1/2
设A(aij)M，定义||A||
|aij|
2
。【7】
i,j1
结合本题数据，利用模型（2）建立规划求解：
优化后的一步转移概率矩阵记为P(2009)P
(20P0*
8)
[114/1510.0.395000.0500
778/819510.]
2111
0.9106
0.8
目标函数：
minfP()
i
P||
*
P||
i1
约束条件：
由MATLAB软件求解得
P*
0.95000.050
0.7889
0.2 ; 111
以2008年基期，预测2009年的概率分布有P(2009)P
(200P*
8)
0.9106150.1P*08(924009)137.494213.5058
则预测2009年发生凝冻的天数为 14天。
6、模型的优缺点分析
本文从整体和局部两个角度分别建模，两个模型的结果大致吻合，说明
Markov链模拟的效果还不错。对于本题，雨凇的出现可以看成是一个随机过程，而影响雨凇出现的因素太多，若直接分析影响因素，会十分麻烦，且由于干预变量会使模型的精度大为降低。
Markov链模型的优点在于它不需要从复杂的预测因子中寻找个因素之间的相互规律，直接通过概率矩阵来预测内部状态的变化。可是Markov链对转移概率矩阵的要求很高，实际中，由于观测数据的限制，误差影响，会造成转移
13
089
概率矩阵精度的降低。因此为保证准确的数据。
Markov链模拟的准确性，需要收集足够多的
参考文献
【1】张波，张景肖.应用随机过程.北京：清华大学出版社，2008年，P74.. 【2】张波，张景肖.应用随机过程.北京：清华大学出版社，2008年，P75.【3】张波，张景肖.应用随机过程.北京：清华大学出版社，2008年，P81.【4】于波, 陈希镇, 华栋.马尔柯夫链在农作物年景预测中的应用.统计与决
策2007年第21期.
【5】历史天气数据http://www.lstqw.cn/lstqw/ 2009年7月23日.【6】于波, 陈希镇, 华栋.马尔柯夫链在农作物年景预测中的应用.统计与决
策2007年第21期.
【7】 [ppt] Chapter1.2 向量范数与矩阵范数—武汉大学精品课程
http://jpkc.whu.edu.cn/jpkcsite/szfx/Data/2009年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=black 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
14
;
procarimadata=text2; identifyvar=freq; run;
2.数据
2002年出现雨凇的天数1月2月3月11月0000000000000000000000000000000000000000000000000000000000000000000000001000100000000000000010001000000000000000
00000000
15
月
0000000000000000000000001111010
122003年出现雨凇的天数1月2月3月11月0000010000000000101010100000000010000000000001000000000000000000000000000000000000000000000000000000000010000000
0000000
0
2004年出现雨凇的天数
1月2月3月11月
000001000
1
0
0
16
12月
0000000000000000000000000000000
12月
000
01001000000000000000000000000000000000000010000010000010010000010010000000
000000
0
2005年出现雨凇的天数月2月3月110000000000000000000000100110
000000000000000000000000000
11月
0000000000
17
0000000000000000000011011111
月
0000110000
11211001001000000000001000000010010010000000000010010000
000001
0
2006年出现雨凇的天数月2月3月0000000000101001001000000000000000000010000000000
1
0
00000000000000000000
11月
00000000000000000
18
000000000000000000000
月
00000000000000000
11201000010010010000000000000000000
000000
0
2007年出现雨凇的天数1月2月3月000000100000000000001000000000000100100100100100100000000000000000000
0000000000000
11月
00000000000000000000000
19
00000000000000
月
00000000000000000000000
1200000000000000
000000
0
2008年出现雨凇的天数月2月3月01001001001001001001001001001001001011010010010010010010010010010010010010011010010
0101
0
000000
11月
000000000000000000000000000000
20
00000000
月
000000000000000000000111000001
112100
21
因篇幅问题不能全部显示，请点此查看更多更全内容