基于马尔科夫链对全国粮食产量的预测--毕业论文.doc

摘要本文主要是运用应用随机过程中的马尔科夫链，通过它的马尔可夫性也就是无后效：通过即要确定过程将来的状态，知道它此时的情况就够了，并不需要对它以往状态的认识。利用马尔可夫链模型,结合实际数据构建二次规划模型,运用Excel与MATLAB统计软件对转移概率矩阵进行了估计和预测,并对从1982年到2012的全国粮食产量进行划分状态，即求出它的转移概率矩阵，然后对以后的年份所处的状态进行预测和分析，也就是N步转移概率矩阵，然后在对2013年的数据进行比对，看是是否正确等关键词：马尔科夫链模型 MATLAB 无后效性 转移概率矩阵 目 录1基本模型与研究思路11.1马尔可夫链模型说明11.2研究思路22模型设计42.1马尔科夫链定义模型42.2 n步转移概率，C-K方程53实际问题的马尔科夫求解63.163.273.384 结论分析10参考文献111基本模型与研究思路1.1马尔可夫链模型说明和随机过程 在当代科学与社会的广阔天地里，人们都可以看到一种叫作随机过程的数学模型：从银河亮度的起伏到星系空间的物质分布、从分子的布朗运动到原子的蜕变过程，从化学反应动力学到电话通讯理论、从谣言的传播到传染病的流行、从市场预测到密码破译，随机过程理论及其应用几乎无所不在。一些特殊的随机过程早已引起注意，例如1907年前后，.马尔可夫研究过一列有特定相依性的随机变量，后人称之为马尔可夫链（见马尔可夫过程）；又如1923年N.维纳给出了布朗运动的数学定义（后人也称数学上的布朗运动为维纳过程），这种过程至今仍是重要的研究对象。虽然如此，随机过程一般理论的研究通常认为开始于30年代。1931年，.柯尔莫哥洛夫发表了概率论的解析方法；三年后，.辛钦发表了平稳过程的相关理论。这两篇重要论文为马尔可夫过程与平稳过程奠定了理论基础。稍后，P.莱维出版了关于布朗运动与可加过程的两本书，其中蕴含着丰富的概率思想。1953年，J.L.杜布的名著随机过程论问世，它系统且严格地叙述了随机过程的基本理论。1951年伊藤清建立了关于布朗运动的随机微分方程的理论(见随机积分)，为研究马尔可夫过程开辟了新的道路；近年来由于鞅论的进展，人们讨论了关于半鞅的随机微分方程；而流形上的随机微分方程的理论，正方兴未艾。60年代，法国学派基于马尔可夫过程和位势理论中的一些思想与结果，在相当大的程度上发展了随机过程的一般理论，包括截口定理与过程的投影理论等，中国学者在平稳过程、马尔可夫过程、鞅论、极限定理、随机微分方程等方面也做出了较好的工作。马尔可夫链可以对动态变化的随机现象进行预测,但应用马尔可夫链进行预测需要满足两个假设,即马尔可夫性。其一是无后效性,就是将来的状态只与现在的状态有关,与过去无关。将来t+1时刻状态仅依赖于第t时刻状态的分布,与过去t-1, t-2,时刻的状态分布及转移状态无关。其二是转移概率稳定性,即马尔可夫理论以固定的转移概率矩阵为基本规律和特征1,要求转移概率矩阵具有相对稳定性2。在长期预测中,转移概率矩阵很难保持不变,因而马尔可夫链预测比较适合短期预测。通过对市场现象的大量观察发现,同类商品的市场占有率分布是一个随时间不断变化的随机过程,并且当期市场占有率只与前一期的市场占有率有关,而与再远期的市场占有率关联甚微3,市场占有率的这一特性与马尔可夫链的无后效性相吻合。经过多年的发展,移动通信市场也基本服从这一规律,尤其从移动客户对运营商选择的心理来看,移动客户选择运营商通常依据其对当前各家运营商的评价,而与再远期的评价关系愈来愈小。另外,马尔可夫的转移概率稳定性包含了对研究期市场和移动客户偏好稳定的特殊要求,而根据移动客户的消费习惯,在运营商的市场政策、政府的监管政策没有特殊变化的情况下,由于移动客户通常存在较高的转网成本,因此维持现有转移概率的可能性较大。本文应用的马尔可夫链是状态集和时间集均为离散的齐次马尔可夫过程。若已知马尔可夫链的一步转移概率矩阵和初始分布,就可得到绝对分布,假设初始市场占有率为S(0)=(s1(0), s2(0), sN(0),转移概率矩阵为P,则m个周期后的市场占有率为：S(m)=S(0)Pm=S(m-1)P(1)1.2研究思路和马尔科夫过程的发展 20世纪50年代以前，研究马尔可夫过程的主要工具是微分方程和半群理论（即分析方法）；1936年前后就开始探讨马尔可夫过程的轨道性质，直到把微分方程和半群理论的分析方法同研究轨道性质的概率方法结合运用，才使这方面的研究工作进一步深化，并形成了对轨道分析必不可少的强马尔可夫性概念。1942年，伊藤清用他创立的随机积分和随机微分方程理论来研究一类特殊而重要的马尔可夫过程扩散过程，开辟了研究马尔可夫过程的又一重要途径。出于扩大极限定理应用范围的目的，马尔科夫在20世纪初开始考虑相依随机变量序列的规律，并从中选出了最重要的一类加以研究。1906年他在大数定律关于相依变量的扩展一文中，第一次提到这种如同锁链般环环相扣的随机变量序列，其中某个变量各以多大的概率取什么值，完全由它前面的一个变量来决定，而与它更前面的那些变量无关。这就是被后人称作马尔科夫链的著名概率模型。也是在这篇论文里，马尔科夫建立了这种链的大数定律。用一个通俗的比喻来形容，一只被切除了大脑的白鼠在若干个洞穴间的蹿动就构成一个马尔科夫链。因为这只白鼠已没有了记忆，瞬间而生的念头决定了它从一个洞穴蹿到另一个洞穴；当其所在位置确定时，它下一步蹿往何处与它以往经过的路径无关。这一模型的哲学意义是十分明显的，用前苏联数学家辛钦（18941959的话来说，就是承认客观世界中有这样一种现象，其未来由现在决定的程度，使得我们关于过去的知识丝毫不影响这种决定性。这种在已知“现在”的条件下，“未来”与“过去”彼此独立的特性就被称为马尔科夫性，具有这种性质的随机过程就叫做马尔科夫过程，其最原始的模型就是马尔科夫链。这既是对荷兰数学家惠更斯（Ch. Huygens, 16291659）提出的无后效原理的概率推广，也是对法国数学家拉普拉斯（P. S. Laplace, 17491827）机械决定论的否定。这里应该指出，尽管拉普拉斯对概率论的早期发展作出过重大贡献，但是他的部分哲学观点是不利于这门学科的深入发展的。十八世纪以来，随着牛顿力学的彻底胜利，一种机械唯物主义的决定论思潮开始在欧洲科学界蔓延，鼓吹最力者就是拉普拉斯。1759年他在巴黎高等师范学院发表了一篇题为概率论的哲学探讨的演讲，淋漓尽致地表达出了这种思想。他说：“假如有人知道了某一时刻支配自然的一切力，以及它的一切组成部分的相对位置，又假如他的智力充分发达，能把这一切数据加以充分的分析，把整个宇宙中从最巨大的天体到最微小的原子的一切运动完全包括在一个公式里面，这样对他就没有什么东西是不确定的了，未来也好，过去也好，他都能纵览无遗。”1812年，拉普拉斯又进一步提出“神圣计算者”的观念，认为这个理想的数学家只须知道世界某一时刻的初始状态，就可以从一个无所不包的微分方程中算出过去和未来的一切状态。换句话说，他认为任意系统在 t t0时的状态 x可由其初始时刻 t0和初始状态 x0唯一决定。这可真是笔判终身、细评流年，数学家可以摆个卦摊了。马尔科夫的概率模型从根本上否定了系统中任一状态 x与其初始状态 x0之间的因果必然性，从而也否定了“神圣计算者”的神话。还应该指出，马尔科夫所建立的概率模型不但具有深刻的哲学意义，而且具有真实的物质背景，在他的工作之前或同时，一些马尔科夫链或更复杂的随机过程的例子已出现在某些人的研究中，只不过这些人没有自觉地认识到这类模型的普遍意义或用精确的数学语言表述出来罢了。例如苏格兰植物学家布朗 ( R. Brown, 17731858) 于1827年发现的悬浮微粒的无规则运动、英格兰遗传学家高尔顿（F.Galton, 18221911) 于1889年提出的家族遗传规律、荷兰物理学家埃伦费斯特 ( P. Ehrenfest, 18801933) 于1907年关于容器中分子扩散的实验，以及传染病感染的人数，谣言的传播，原子核中自由电子的跃迁，人口增长的过程等等，都可用马尔科夫链或过程来描述。也正是在统计物理、量子力学、遗传学以及社会科学的若干新课题、新事实面前，决定论的方法显得百孔千疮、踵决肘见。有趣的是，马尔科夫本人没有提到他的概率模型在物理世界的应用，但是他利用了语言文学方面的材料来说明链的性质。在概率演算第四版中，他统计了长诗叶甫盖尼奥涅金中元音字母和辅音字母交替变化的规律：这是长诗开头的两句，意为：“我不想取悦骄狂的人生，只希望博得朋友的欣赏。”诗人那火一般的诗篇在数学家那里变成了一条冷冰冰的锁链：在这条锁链上只有两种链环，C代表辅音、 代表元音（为了使问题简化起见，不仿把两个无音字母算作辅音）。马尔科夫分别统计了在C后面出现C和 的概率p和1p，以及在 后出现C和 的概率q和1q，把结果与按照俄语拼音规则计算出的结果进行比较，证实了语言文字中随机的（从概率的意义上讲）字母序列符合他所建立的概率模型。完成了关于链的大数定律的证明之后，马尔科夫又开始在一系列论文中研究链的中心极限定理。1907年他在一种不平常的相依试验中证明了齐次马尔科夫链的渐近正态性；1908年在一个链中变量和的概率计算的极限定理推广中作了进一步的推广；1910年他发表了重要的论文成连锁的试验，在其中证明了两种情况的非齐次马尔科夫链的中心极限定理。与此同时他在一些假定的前提下证明了模型的各态历经性，成为在统计物理中具有重要作用的遍历理论中第一个被严格证明的结果。遍历理论亦称ergodic理论, 是奥地利物理学家玻耳兹曼（L. Boltzmann, 18441906) 于1781年提出来的，其大意是：一个系统必将经过或已经经过其总能量与当时状态相同的另外的任何状态。马尔科夫链的引入，在物理、化学、天文、生物、经济、军事等科学领域都产生了连锁性的反应，很快地涌现出一系列新的课题、新的理论和新的学科，并揭开了概率论中一个重要分支随机过程理论蓬勃发展的序幕如果可以确定粮食来年状态的所有可能状态及其初始分布,进而确定变量的状态转移概率4,建立马尔可夫预测模型,根据模型求出绝对分布,就可预测短期的粮食状态。其中转移概率矩阵是马尔可夫链理论预测的核心问题5,但实际情况中,并不总是能从市场调查直接得到转移概率矩阵。本文采用二次规划模型,结合实际数据求解转移概率矩阵。过去,该方法计算过程较为复杂,有时手工计算难以完成,难以将其推广应用;也有些学者利用LINGO软件进行求解,但通用性不强,针对每个问题都必须重新编程,对于没有计算机编程基础的预测者而言并不实用6。因此本文应用内置于Office产品Excel中的规划求解工具实现马尔可夫转移矩阵计算,不仅操作简便通用,而且可大大提高预测工作的效率和准确度。2.1马尔科夫链定义模型随机过程称为马尔科夫链，若它只可以取有限或可列个值，并对任意状态有： 在上式中 表示过程 处状态 ，称为0,1,2，为该过程的状态空间，记为S。 当马尔科夫链的转移概率 只与 有关，而与n无关时，称为齐次马尔科夫链，否则称为非齐次。其中pij为状态i到状态j的转移概率。我们将 排成一个矩阵的形式，令P= 称P为转移概率矩阵，一般简称为转移矩阵，由于概率非负，且过程必须从一过程转移到另一过程，所以有： （1） （2） 2.2 n步转移概率，C-K方程 称条件概率 为马尔科夫链的n步转移概率，相应地称 为n步转移概率矩阵。对一切 有（1） （2）3.实际问题的马尔科夫预测 例如，考虑中国几十年的粮食产量分为三个状态，即“丰收”、“平收”和“欠收”。记E1为“丰收”状态，E2为“平收”状态，E3为“欠收”状态。下图给出了全国19822012年期间全国粮食产量的状态变化情况，1981年的全国粮食总产量为32502（单位：万吨） 年份 粮食产量 增长速度 状态 1982年：35450 8.32%丰 E11983年：38728 8.46%丰 E11984年：40731 4.92% E21985年：37911-7.44%欠 E31986年：39151 3.17% E21987年：40473 3.27% E21988年：39404-2.71% E2 1989年：40755 3.31% E21990年：44624 8.67%丰 E11991年：43529-2.52% E2 1992年：44266 1.66% E2 1993年：45649 3.03% E21994年：44510-2.56% E21995年：46662 4.61% E21996年：50454 7.52%丰 E11997年：49417-2.10% E21998年：51230 3.54% E21999年：50839-0.77% E22000年：46218 -10.00%重欠 E32001年：45262-2.11% E22002年：45711 0.98% E22003年：43067-6.14%欠 E32004年：46947 8.26%丰 E12005年：48401 3.00% E22006年：49746 2.70% E22007年：50150 0.81% E22008年：52850 5.11% E12009年：53082 0.44% E22010年：54641 2.85% E22011年：57121 4.34% E22012年：58957 3.11% E23.1首先计算状态转移概率矩阵假定某一个事件的发展过程有n个可能的状态，即E1，E2， ，En。记为从状态转变为状态的状态转移概率，则矩阵 从表3.7.1中可以知道，在6个从E1出发（转移出去）的状态中，有1个是从E1转移到E1，有5个是从E1转移到E2，有0个是从E1转移到E3，所以按照上述同样的办法计算可以得到所以，中国粮食产量的状态转移概率矩阵为 E2E1 0.167 0.136 0.728 0.833 0 0.136 0.667E3 0.333 03.2进行预测计算状态概率表示事件在初始（k0）状态为已知的条件下，经过k次状态转移后，在第k 个时刻（时期）处于状态的概率。根据概率的性质，显然有 从初始状态开始，经过k次状态转移后到达状态这一状态转移过程，可以看作是首先经过（k-1）次状态转移后到达状态，然后再由经过一次状态转移到达状态。根据马尔可夫过程的无后效性及Bayes条件概率公式，有： 若记行向量，则由（3.7.7）式可以得到逐次计算状态概率的递推公式： 式中：为初始状态概率向量。3.3利用Matlab编程进行马尔可夫预测计算以1982年的农业收成状态为初始状态，预测今后年（即20132018） p=0.167 0.833 0;0.136 0.728 0.136;0.333 0.667 0;for i=1:1:6 y=piend y = 0.1670 0.8330 0 0.1360 0.7280 0.1360 0.3330 0.6670 0y = 0.1412 0.7455 0.1133 0.1670 0.7340 0.0990 0.1463 0.7630 0.0907y = 0.1627 0.7359 0.1014 0.1607 0.7395 0.0998 0.1584 0.7378 0.1038y = 0.1610 0.7389 0.1001 0.1606 0.7388 0.1006 0.1614 0.7383 0.1003y = 0.1607 0.7388 0.1005 0.1608 0.7387 0.1005 0.1608 0.7388 0.1004y = 0.1608 0.7387 0.1005 0.1608 0.7388 0.1005 0.1608 0.7388 0.1005 p=0.167 0.833 0;0.136 0.728 0.136;0.333 0.667 0; x=1 0 0; for i=1:1:6 y=x*piendy = 0.1670 0.8330 0y = 0.1412 0.7455 0.1133 y = 0.1627