全国数学建模大赛C题

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

2、置报名号的话）：139C01所属学校（请填写完整的全名）：浙江工贸职业技术学院参赛队员（打印并签名）：1.郑济明2 .王庆松3 .朱松祥指导教师或指导教师组负责人（打印并签名）：王积建日期：2012年9月幽日赛区评阅编号（由赛区组委会评阅前进行编号）：编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）： 赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）： 全国评阅编号（由全国组委会评阅前进行编号）：脑卒中发病环境因素分析及干预摘要关键词：一、问题重述21世纪人类倡导人与自然和谐发展，环境因素成为影响健康的重要因素。脑卒中 （俗称脑中风）就是与环境因

3、素紧密相关且威胁人类生命的疾病之一。这种疾病的诱 发已经被证实与环境因素有关，其中与气温和湿度存在着密切的关系。对脑卒中的发 病的环境因素进行分析，其目的是为了进行疾病的风险评估，对脑卒中高危人群能够 及时采取干预措施，也让尚未得病的健康人，或者亚健康人了解自己得脑卒中风险程 度，进行自我保护。同时，通过数据模型的建立，掌握疾病发病率的规律，对于卫生 行政部门和医疗机构合理调配医务力量、改善就诊治疗环境、配置床位和医疗药物等 都具有实际的指导意义。现从中国某城市各家医院2007年1月至2010年12月的脑卒中发病病例信息以及 相应期间当地的逐日气象资料（Appendix-C2）和数据（见 Ap

4、pendix-C1 ）。需解决一 下几个问题：问题一：根据病人基本信息，对发病人群进行统计描述。问题二：建立数学模型研究脑卒中发病率与气温、气压、相对湿度间的关系。问题二：查阅和搜集文献中有关脑卒中高危人群的重要特征和关键指标，结合 1、 2中所得结论，对高危人群提出预警和干预的建议方案。二、问题分析脑卒中（俗称脑中风）作为威胁人类生命的疾病之一，并且病发的人群受环境因 素的影响不断扩展。对脑卒中人群及受环境因素的影响分析来对疾病的风险评估，对 脑卒中高危人群能够及时采取干预措施成为一项无疑是一项十分复杂的系统工程。对于问题一，利用中国某城市各家医院2007年1月至2010年12月的脑卒中发病

5、 病例信息以及相应期间当地的逐日气象资料 （Appendix-C2）和数据（见Appendix-C1 ）。 通过excel对已知数据进行统计整理，再利J用matlab程序对脑卒中病发者的性别指数、 年龄指数、职业指数、月份指数进行合理的统计得出相应数据比率。三、模型假设3.1 模型假设：1）发病病例的信息中，若两个病例的信息相同，则视为不同的两个人；2）以诊断报告时间为准来统计发病人群的数量；3）导致脑卒中发病的内在原因只与性别、年龄、职业有关；4）导致脑卒中发病的外在原因只与气压、温度和湿度有关；5）气压、温度和湿度之间具有相关关系；6）月平均气压、月平均最高气压、月平均最低气压具有相关关系

6、；7）月平均温度、月平均最高温度、月平均最低温度具有相关关系；8）月平均湿度、月平均最高湿度、月平均最低湿度具有相关关系；9）关于环境因素如气压、温度和湿度的观测数据都是准确可靠的；10）按照国际惯例1,发病率以10万人群的发病人数来表示。但由于本题是研究 某地区的发病人数，并没有与其它地区比较，所以在本题分析中，发病率以发病人数 来表示。四、符号说明定义1,月平均气压是日平均气压的平均值。月平均最高气压是日平均最高气压 的平均值。月平均最高气压是日平均最高气压的平均值。定义2,月平均温度是日平均温度的平均值。月平均最高温度是日平均最高温度 的平均值。月平均最高温度是日平均最高温度的平均值。定

7、义3,月平均湿度是日平均湿度的平均值。月平均最低湿度是日平均最低气压 的平均值。Ni表示男性病例总数，心表示女性病例总数，N表示总病例数；五、模型的建立及求解5.0发病人群数据的预处理根据已知题意给出的中国某城市各家医院 2007年1月至2010年12月的脑卒中发病病 例信息以及相应期间当地逐日气象资料，进行如下数据预处理：1）以2007年1月至2010年12月的一共48个月的脑卒中发病病例为准，其他时间数 据应当删除，一共得到58925个病例.2）如果病例的信息中，年龄与职业不符（例如：12周岁是老师）、诊断时间不详、数据明显出错的都不应该考虑在统计范围之内，应当删掉。3）3.1 对发病人群

8、的统计分析（问题1）3.1.1 性别分析1）性别差异性简单分析男、女性病发比例为N4=1=1,2 （1）N其中，i=1表示男性，i=2表示女性。经统计，N =58925, Nn =31832, N12 = 27093 ,代入（1）得男、女病发比例分别为 54.02%45.98%（matlab程序见附录1）。可见男性在脑卒中的病发者要大于女性脑卒中病发人 数。2）单因素方差分析2逐月统计男女病例人数，考察在相同时间点上男女人群发病人数是否有显着差异，给定显着性水平a =0.05,分析结果为F =5.54,对应的p = 0.0206 < 0.05 （见图1）,又查表得 F2（r-1,n-r）

9、 = F）.05 （2 - 1,48- 2）=5（1,46米 4.0施于 F > Fo.o5（1,46）,所以脑卒中发病男女人群有显着差异（matlab程序见附录2）。图1男女发病人群的单因素方差分析结果1）简单分析不同年龄阶段发病比例为N2i-X2i = i =123,4,5,6 （2） N其中，i =1,2,3,4,5,6 分别表示 “40 岁以下”、“40-50”、“50-60”、“60-70”、“70-80”、“80以上”。经统计，将不同年龄阶段脑卒中病发者人数代入（2）式，得到不同年龄 阶段脑卒中病发者比例，见图2。（matlab程序见附录3）图2不同年龄段发病人群比例图由图2

10、可以看出在50岁以下的人口中脑卒中病发的人数比例较小，70-80之间脑 卒中的比例最为严重，80岁以上的人脑卒中较为严重，所以高龄的人是发生脑卒中的 高危人群，我们应当高度关注。2）单因素方差分析根据图2结果，剔除“ 40岁以下”和“ 40-50”年龄段，对其余4个年龄段进行 单因素方差分析，逐月统计不同年龄段发病人群人数，考察在相同时间点上不同年龄 段发病人群人数是否有显着差异， 给定显着性水平u =0.05,分析结果为F =45.6,对应的p =0.0000 <0.05 （见图3）,所以脑卒中不同年龄段发病人群有显着差异 （matlab程序见附录 4 （tongji6.m 和 ton

11、gji60.m ）。图3不同年龄段发病人群的单因素方差分析结果1）简单分析不同职业发病比例为N3i . 一 -X3i =, i = 1,2,.,9 （3）N其中，i =123,4,5,6,7,8,9分别表示“农民”、“工人”、“退休人员”、“教师”、“渔民”、“医务人员”、“职工”、“离退人员”、“其它职业”。经统计，不同职业脑卒中病发 者的比例，见图4。（matlab程序见附录4）图4不同职业病发者比例由图4得出农民、工人、退休人员、其他职业的人员患脑卒中的比例偏高，说明 了职业也是患脑卒中的重要因素。2）单因素方差分析根据图4结果，对农民、工人、退休人员进行单因素方差分析，给定显着性水平a

12、 =0.05,分析结果为F =95.36,对应的p = 0.0000 < 0.05 （见图5）,所以脑卒中不同年龄段发病人群有显着差异（matlab程序见附录5 （tongji7.m 和tongji70.m ）。图5不同年龄段发病人群的单因素方差分析结果1）简单分析定义季节指数为x4i = -411=1,2,.,12 （4）S其中，S4i为第i月的平均人数，S为48个月的月平均人数。经统计，不同月份脑 卒中病发者的比例，见图6。（ matlab程序见附录6）图62007年-2010年各月季节指数由图6看出在五、六、七月份为脑卒中高发期，一、二月为低发期 。5.2发病率与气压、气温、相对湿

13、度间的关系分析（问题2）由于题目提供了环境因素（气压、温度和湿度）的 8个变量，根据假设5）8）, 这8个变量间具有明显的显着相关关系，所以必须做降维处理，把8个变量整合成立不相关的少数几个变量，然后再寻找发病率与这少数几个变量的关系式。这需要进行 主成分分析。 31）基本原理主成分分析是把多个变量转化为少数几个新综合变量的一种多元统计方法，其基 本思想就是在保留原始变量尽可能多的信息的前提下达到降维的目的，从而简化问题 的复杂性并抓住问题的主要矛盾.其手段是将原来众多的具有一定相关性的变量重新 组合成新的少数几个相互无关的综合变量（也叫抽象变量），来代替原来变量，这些新 的综合变量称之为主成

14、分.一般地说，利用主成分分析得到的主成分与原来的变量之间有如下基本关系：（1）每一个主成分都是各原始变量的线性组合.（2）主成分的数目大大少于原始变量的数 目.（3）主成分保留了原始变量的绝大多数信息.（4）主成分之间互不相关.据此我们 建立数学模型.2）数学模型在一个统计问题中，假设我们收集到 n个样品，每个样品观测到p个变量（记为Xi,X2,|xp,为简单起见，可以设Xi均值为0,方差为1, （1 <i < p ）,构成一个nM p阶的样本原始资料阵X = %.ij n：p主成分分析的目的在于利用p个原始变量（X1,X2,m,xp）构造少数几个新的综合变量，使得新变量为原始变量

15、的线性组合，新变量互不相关，新变量包含p个原始变量的绝大部分信息.这样定义X1,X2,| ,Xp为原始变量，y1,y2,IM , ym（m w p）为新的综合变量指标，每一个新综合变量指标是 p个原始变量的线性组合：yi =aiiXi aI2X2 . IH . aipXpY2 =a2iXa22X2a2Pxp /l(5)ym amiXi am2X2 a ampXp同时要求满足以下几个条件：(1) yi与“相互无关；(2) yi是Xi,X2，川，Xp的一切线性组合中方差最大者；y2是yi与不相关的Xi, X2，川，Xp的所有线性组合中方差最大者;|,ym是Zi,Z2|,Zm分别都不相关的Xi, X

16、2，川，Xp的所有线性组合中方差最大者.则新、第二、第 m主成分.变量yi, 丫2,川，ym分别称为原变量Xi, X2,| , Xp的第从以上的分析可以看出，主成分分析的实质就是确定原来变量Xj( j = i,2,| , p)在诸主成分yi(i =i,2/|,m)上的系数aij(i=i2|,m; j=i,2JH,p).从数学上可以证明，他 们分别是p个原始变量(Xi,X2,|,Xp)相关矩阵的前m个具有较大特征值所对应的特 征向量，而各个新综合变量 X的方差var(yi)恰好是相应的特征值 片.各主成分的方差 贡献大小按特征根顺序排列，是依次递减的，即九i之九2川II之九p之0 .其几何意义是

17、：主成分分析相当于对原坐标轴做一次旋转变换，使得新坐标系的第i轴对应于数据变易的最大方向，第2轴与第i轴正交，且对应于数据变易的第二大方向，依次类推.3)基本步骤(D确定分析变量，收集原始数据；设原始数据矩阵为 X =(Xj )n即其中Xj表示第i个样品(对象)在第j个变量上的取值。(2)在进行主成分分析之前，要检验该样本矩阵是否适合于主成分分析.KMO检验是检验变量之间偏相关关系的统计量，用于检验变量间的偏相关系数是否过小.KMO统计量越接近于i,说明各变量间的偏相关系数越大，KM佻计量大于0.9 ,效果最好； 如果统计量小于0.6 ,则不适合于做主成分分析.Bartlett 球形检验是检验

18、相关矩阵是 否是单位矩阵，即各变量是否各自独立.(3)对原始数据进行标准化，即令(6)* Xij - XjXij 二sj其中Xj,Sj分别为第j列元素的样本均值和样本标准差，即 则X* =(X*)n冲为标准化的样本资料库(4)由标准化后的数据矩阵求协方差矩阵Z ,或者由原始数据矩阵求相关系数矩阵R.这两种方法结果相等.本文采用直接计算原始数据的相关矩阵的方法(对于数量 级差别较大或者有量纲的数据宜适用).设原始数据X的相关系数矩阵为Tir12 . 小、RM r2 1r22 r2p.Jnlrn2 . rnp ,rj (i, j =1,2,|, p)为原变量为与为的相关系数，其计算公式为rrj(8

19、).一 (Xki - X )(Xkj Xj) k=1nn- (xki - x).=(xkj - Xj) k=4k=4(5)计算R的特征根和特征向量;根据特征方程RE-R =0得R的特征根为九i(i =1,2,，p),将特征根按照从大到小的顺序排列，排列后的特征根不妨仍然表示为 A之九2'与至0.同时可得对应的特征向量U1 ,U2,.,Up,将他们标准正交化，5 ,U2,.,Up称为主轴.(6)计算所有变量的方差贡献率及累计方差贡献率；九i的方差贡献率为e=4M=1,2川,p) (9)- 1i 1A的累计方差贡献率为"1kEi =Jsf(i =1,2,|, p, k = 1,2

20、,.,m,m< p) (10)i W(7)确定主成分的数目m.方法有：一般取累计贡献率达 85395%勺主成分；选用所有人土1的主成分；累计特征值乘积大于1的主成分；画出特征值变化曲线，以转折点位置为标准判断.本文采用累计贡献率达85395%勺主成分.(8)确定主成分函数表达式模型.设m个主成分对应的特征向量分别为T一A、&、.、Am,其中Aj=(aja2j.apj) ,0<j表示Aj的第k行的元素，则第j个主成分yj的函数表达式为T X2yj =(A )=(即 a2j.IxpJX2papj )=£ akjXk (11).k1lxp>(9)提炼主成分yj的抽

21、象意义.由Xk与yj的相关系数b. j的大小可以确定yj主要与哪几个变量显着相关，然后根据这几个变量的实际意义提炼yj的抽象意义.(10)检验主成分模型.根据n个样本的m个主成分的函数值，通过计算m个主成 分y1, y2,., ym的相关系数就可以检验 m个主成分是否线性无关.如果两个主成分的相 关系数为0,则说明这两个主成分线性无关，模型有效；否则线性相关，模型无效.(11)求主成分函数值。将各样本标准化数据Xk代入(7),可以求得各样本的第j个主成分yj的函数值.4)模型求解(1)收集原始数据矩阵X.本文选取了某地区的月平均气压的平均值、月最高气压的平均值、月最低气压的 平均值、月平均气温

22、的平均值、月最高气温的平均值、月平均气压的平均值8项指标,并分别记为X1, X2X8.每个指标有48个数据(见附件1)。使用SPSSa件进行求解(见附录7)。(2)将原始数据标准化，(SPS汕部计算).(3)求原始数据的相关系数矩阵 R ,如图7所Correlation MatrixXiv2*3x4x5暄x7时Cflrnelalio(i1000998999-94Qr3g 376 303x2.9901.000993-.951-S5Q-J950. 3B5-.313x3g邨9931 WO-921921-.917加。265x4-.951-.9211 OOO.aea弼424必-93995。-92199S

23、1 QOD991355263x6由50917MB的11.000,461SBOk7-370-385-3604243854611 000822喟-303-313286,323263JOO.0221.000示.图7相关系数矩阵图8因子分析检验图从图8看出，表格的第一行为检验变量间偏相关程度的KMOS计量，其值在0.6之上才适合做主成分分析，效果显着，如果小于0.6,效果不显着，不适合做主成分分析。下面的三行为球形检验的结果，球形检验原假设的变量是不相关的，显然只有 拒绝原假设的情况下数据才适合做因子分析。本例中KMOB为0.720,球形检验显着， 两个条件都满足，变量间相关程度大，适合做因子分析。(

24、4)计算矩阵R的特征根、各因子的方差贡献率及累计方差贡献率，并确定主成 分的个数.如图9所示。图9R特征值及其累计方差贡献率从图9中可以看出，第一、第二主成分对方差的累计贡献率达到95.461%,它们分别对应着原样本数据点数据变异的最大、次大方向，是原变量系统的一个最佳整合， 从而我们可以以95.461%的精度将变量白有效维数从8维降至2维.因此可以将前2个 因子作为主因子.（5）确定主成分函数表达式模型，因子得分系数矩阵如图10所示.图10因子得分系数矩阵设2个主成分分别为y1，y2,则建立模型为乂 =-0.184x10.183x2 0.185&+0.175x4 +0.186x5 +

25、0.164x6 0.096x7 0.121x812）y =0.048为 +0.042x2 +0.057x3 -0.020x4 -0.055x5 +0.015x6+0.536x7 + 0.567x8其中x1，x2,x8.均为原变量经过均值为0,方差为1标准化后的变量.（6）对主成分yhy2的意义进行解释。图11给出了原变量与第1、第2主成分的 相关系数.图11旋转后的因子载荷矩阵第一主成分力 ,与原变量不（平均气压的平均值）、x （最高气压白平均值）x3 （最低气压的平均值）、% （最低气温的平均值）的相关系数的绝对值都超过了0.948,因此它是一个反映气温和气压的综合因子，我们称之为气压温度因

26、子.第二主成分y2,与原变量x （月平均相对湿度的平均值）的相关系数为0.925、x8 （月最低相对湿度的平均值）的相关系数为0.948,其余的都不超过0.266,因此它是 一个反映相对湿度的因子，称为湿度因子.（7）计算2个主成分的函数值.将48个经过标准化的数据x*代入模型yj，可以得到48个地区的主成分yj的函数 值，结果如表1所示.表148个月对应的主成分函数值在舁 厅Py1y2z在舁 厅py1y2z1-1.68170.548548725-1.601-0.0619462-0.7683-0.208252826-1.1211.777993-0.4507-0.216678427-0.6750

27、.00767394r 0.152-1.5909858280.1144-0.9178615r 1.1598-1.77581082290.9219-2.09885560.92590.75541137301.34170.05847887I 1.52390.20971043311.35210.471575281.5284-0.42721080320.94021.664496290.56751.00551206330.35751.511582810-0.13180.1921265340.2595-0.72776911-0.6524-0.9464120235-1.0911.166766712-1.1893

28、0.5304129736-1.3310.315362313-1.5177-0.4902120937-1.4530.257784114-1.3599-1.3556171638-1.0830.8107128015-0.1654-2.0187178739-0.82301845-1.03159840-0.401-0.076155317:0.9955-1.31491755410.7358-0.6591715180.67641.86591533420.78240.70991618191.6182-0.25252078431.24550.6445162120P 1.14590.73

29、561322441.3537-0.1331716210.5359158160.9467164722-0.18990.9149133146-0.2460.4312170923-0.90770.5766147547-0.531-0.633161524-1.036-0.9904147248-0.593-1.4721719(8)检验主成分模型由于主成分分析的4个条件中的前3个(每一个主成分都是各原始变量的线性组 合；主成分的数目大大少于原始变量的数目；主成分保留了原始变量的绝大多数信息)， 只要检验4个主成分是否相关即可.由步骤(7)计算的2个主成分的得分矩阵Y=(y0 )4

30、叼，求矩阵Y的协方差矩阵如 图12所小.图12因子得分的协方差矩阵从图12可以看出，主成分得分的协方差矩阵为单位矩阵，说明提取的2个主成分是互不相关的.满足假设的条件，模型和结果有效.5.2.2多元非线性回归分析4将48个月的发病率作为因变量，记作z,将发病率的48个数据填入表1中。下面寻找发病率z与主成分。2的关系式，这需要使用多元非线性回归分析方法。经过反复试验探索，找到的非线性回归模型为z = expy十 a2y2 + a3' 十a4y2 + a5十 % I (13)y y2 y1(1 y2) y1其中，a1 =0.1281,a2 = -0.0472,a3 = 0.0273 ,

31、a4 = 0.0021,a5 =0.0407, a6 = 7.0361。模型检验的p =0.0399 <0.05 ,说明模型有效。平均绝对相对误差为26.38%。5.2.3结果分析1)从非线性回归模型(12)可以得到以下结论：(1)由4>0可得，发病率与气压温度因子具有正相关性；(2)由a2<0可得，发病率与湿度因子具有负相关性；(3)由a1 >同可得，气压温度因子比湿度因子对于发病率的影响显着；2)从主成分模型(11)可以得到以下结论：由第1个方程可知：(1)由XhX2,X3的系数为负值可得，气压温度因子与月平均气压、月平均最高气压、月平均最低气压具有负相关性；(2)

32、由X4,%,X6的系数为正值可得，气压温度因子与月平均温度、月平均最高温度、月平均最低温度具有正相关性；(3)由X7,X8的系数为负值可得，气压温度因子与月平均相对湿度、月平均最低相对湿度具有负相关性；由第2个方程可知：(4)由XhX2,X3的系数为正值可得，湿度因子与月平均气压、月平均最高气压、月平均最低气压具有正相关性；(5)由X4<0,% <0,X6>0可得，湿度因子与月平均温度、月平均最高温度成负相关性，与月平均最低温度具有正相关性；(6)由X7,X8的系数为正值可得，湿度因子与月平均相对湿度、月平均最低相对湿度具有正相关性；5.3高危人群预警分析(问题3)首先预测2

33、011年的气象状况，然后预测未来2011年的高危人群的发病率和发病 时间，最后提出预警和干预措施。1)气压、温度状况预测观察月平均气压X月平均最高气压X2、月平均最低气压X3、月平均温度X4、月平均最高温度X5,月平均最低温度X6的历史数据，发现随时间做周期性变化，于是建立余弦函数模型X = Acos(®t+e)+B ,利用过去36个月的历史数据进行参数估计， 然 后使用2009年的12个数据进行预测，评估误差并检验模型的可靠性。最后预测出2010 年的12个月的数据。建模结果见 表2. ( MATLA程序见附录8)。2011年12个月的预 测结果见表3.表2模型参数估计结果变量表i

34、A式参数拟合值平均相 对误差0.15%0.14% 0.15%11.77%7.75% ;34.71%表32011年12个月的气压、温度预测值月份123456平均温度4.99785.29698.666214.211120.459625.7531平均压力1026.11024.31020.11014.71009.61006最高压力1018.21023.71027.710291027.41023.2最低压力1014.71019.51022.61023.21021.11016.9最高温度23.742617.530212.2499.2949.445612.6638最低温度16.02659.69374.389

35、21.52791.87315.3328月份789101112平均温度28.686228.480325.190119.68913.43748.0948平均压力1004.91006.71010.91016.21021.41025.1最高压力1017.61012.11008.11006.61008.21012.3最低压力1011.71006.91003.81003.21005.21009.4最高温度18.098524.314429.669732.7532.741729.647最低温度10.983917.319122.648425.550225.250321.82892)湿度状况预测观察月平均湿度X7

36、、月平均最低湿度的历史数据X8,发现它们随时间做平稳性波 动，于是建立马尔克夫模型，利用过去 36个月的历史数据进行建模，然后预测出第37个数据。采用“新陈彳t谢”思想，把第 37个数据加入建模序列，并同时去掉最老 的第1个数据，保持数据“等维”，建模并预测出第38个数据，如此滚动预测，直至 预测出2009年的12个数据，并做误差分析，检验模型的可靠性。最后预测出2011年的12个数据。(1)自相关系数原始序列X(0) =X(1),X(0) (2)，,X(0) (n)的各阶自相关系数反映已知数据对未来数据的影响程度.各阶自相关系数为n -w(13)、x(0)(k)-X(0) X(0)(k w)

37、-X(0) km 一一n2X X(0)(k)-X(0)2k=1 -式中,nX(0)=-E X(0)(k) (14) n km对各阶自相关系数归一化得,- r”0w = t, w = 1,2，,t (15)rw 1为可作为各阶步长的马尔柯夫链权重,t是按预测需要计算的最大阶数，一般取rw主0.3.根据rw 20.3可以确定转移步数w.(2)加权马尔柯夫模型状态划分。设划分的m个湿度区间为其中，即尽可能小，a2m尽可能大.，如果则表明第k年的相对误差处于第i种状态.状态转移概率矩阵的构造。设 w步转移概率为pj(w),记:m(w)端工方，小12.小(15)其中，m：w)表示状态Ei经过w步转移到状

38、态Ej的次数，Mi为状态Ei出现的次数.由于数据序列最后的状态转向不确定，故计数Mi时要去掉数据序列中最末的 w个数据(也就是只考虑前面的n -w个数据).由p(w)构成的矩阵称为w步转移概率矩阵，记作(16)已知每一步的概率转移矩阵和每一步的初始状态，则马尔柯夫链就可以确定 预测值计算选取距离预测年最近的t(tEm)个年份，按照距离预测年由近到远，转移步数w分别为1,2,t ,以这几年的相对误差所对应的状态为初始状态，不妨设第1,2,t年所对应的初始状态分别为Ell E22).,寸)其中，w'wl,2,m.例如，当2'=5时，E2(2)=E 5(2)，说明距离预测年第2年的状

39、态是第5状态.在转移步数w对应的转移矩阵从而组成新的概R(w)中,取起始状态eWw)所对应的行向量pWMpWZ, pwX., pwX),率矩阵''湍)p221)第P22P(1m、 p22m(17)Pt(t1)Pt('t2)pttm ，将矩阵R加权得记;)询 p(121Mm(18)将矩阵R0按列求和得tw twtp= R,P2,，Pm =VwpWW）J 犷肾,，、 w 4w 1w 1找出向量p的最大分量得Pm = maxiR, P2,Pm）,M il,2,分量Pm所对应的状态Em就是预测年的状态，则该年度的预测值为al M.,m?(19)(20)(21)（3）计算过程和结

40、果以预测2010年第1月的数据为例。利用2007、2008、2009年的36个月的历史数 据进行建模。自相关系数。以根据（13）、（14）、（15）式计算可得各阶的自相关系数，确定最大滞后阶数w = 2 .各阶自相关系数及权重见 表4.表4自相关系数及权重120.47370.24140.66250.3375划分的6种状态区间，见表4.表4各个状态区间状态编号状态区间0,60)60,65)65,70 )70,75)75,80 )80,100构造转移概率矩阵如果有的状态不能从统计表中得到转移概率，则假定它未来转移到各个状态的概率都相等，即都等于 -.根据（15）可得1步和2步内的转移概率矩阵分别见

41、 表5和 m表6.表51步转移概率矩阵R2/31/3000001/52/51/51/5011/81/41/803/81/801/73/71/72/70001/94/91/31/91001/31/301/3表62步转移概率矩阵R1/31/31/30001/5001/52/51/511/81/805/81/8002/74/702/71/701/81/201/41/81001/302/30组成预测年份的新转移概率矩阵.选择离预测年最近的2个年份，转移步数分别为w=1,2,根据(17)式得预测年的转移概率矩阵R37 ,见表7.表7月 份状 态步 长权重预测年的转移概率矩阵R37概率来源12345636

42、210.66250r 1/73/7 12/71/7r 035320.337501/81/201/41/8加权求和00.136 80.45270.09460.273 70.042 2确定预测年份的状态.预测年的状态向量的最大分量值为 0.4527 ,对应的状态 为第3状态，即第37个月的湿度将处于第3状态，湿度67.5.将第37个月的湿度值67.5放入序列中，同时去掉第1个月的湿度数据，重新构 建马尔柯夫链，得第38个月的湿度。以此类推，可得2010年12个月的湿度值，见表 8.平均绝对相对误差为7.97%,可靠性高.表82010年湿度预测月份123456789101112实际x767.577.

43、572.567.577.572.567.577.572.567.577.572. 5模x770.974.670.869.568.376.475.571.577.274.268.962.2481798893绝对相对误4.853.782.392.9813.35.1010.68.276.099.1412.416.5差%2%9%8%0%表82011年湿度预测月 份12345678910111272.77.72.67.77.72.67.77.72.67.77.72.55555555555557.57.57.57.57.57.57.57.57.57.57.57.555555555555将预测得到的2011

44、年12个月的气象数据代入主成分模型(12)和多元非线性回归模型(13),计算得到2011年发病率的发展趋势，如 图13所示。图132011年发病率趋势从图13可知，2011年发病率的发展趋势是，15月保持平稳态势，从6月份开 始呈现上升态势，至8、9、10月份到达高发期，之后，11和12月份下降进入低发期。查阅文献5-9,发现以下信息：1)脑卒中，又称中风或脑血管意外，它包括脑出血、蛛网膜下腔出血、脑梗死和短 暂性脑缺血发作等急性脑血管病，是一组突然起病，以出现意识障碍和局灶性神经功能 缺失为共同特征的急性脑血管病。2)它具有以下特征:(1)发病率高;(2)致残率高;(3)死亡率高;(4)复发

45、率高3)它的关键指标为：(1)头晕。(2)肢体麻木。(3)暂时性吐字不清或讲话不灵。 (4)肢体无力或活动不灵。(5)与平时不同的头痛。(6)不明原因突然跌倒或晕倒。 (7)短暂意识丧失或个性和智力的突然变化。(8)全身明显乏力，肢体软弱无力。(9) 恶心呕吐或血压波动。脑卒中病人大部分是“三高”、有家族病史与前科老年人。建议措施：(1)卫生机构在2011年的秋季前做好迎接高发病人的准备。(2)预防 为主，多锻炼，多吃蔬菜，不吸烟，不喝酒，防“三高”。六、需要进一步研究的问题对高危人群如男性人群、60岁以上人群、农民群体进行预测，做好预警和干预。七、模型评价7.1 模型优点：1)主要因素并进行

46、了定量分析，7.2 模型缺点：1)当考虑定性相关因素较多时，八、参考文献1周晓平，杨进.脑卒中发生时气节规律及其气象医学原理探讨.中医杂志，2006年6 月，第47卷第6期.2谢文龙，尚涛.SPSS统计分析与数据挖掘.北京：北京电子工业出版社，2012.13韩中庚.数学建模方法与应用.北京：高等教育出版社，2005.6(2006重印).4陈在余，陶应虎.统计学原理与实务.北京：清华大学出版社，2009.3复发高危人群，低温危害高温危害九、附录附录1.统计发病人群的性别比例建立M-文件：流计发病人群的性别比例A=xlsread('c 题数据 1.xls','sheet1&

47、#39;,'A2:F58926');%读入发病人群信息，58925 行-6列a=size(A);B=zeros(1,2);fori=1:a(1)ifA(i,1)=1B(1,1)=B(1,1)+1;% 男性病人的人数elseB(1,2)=B(1,2)+1;% 女性病人的人数endendB=B/a(1),%性别比例附录2.方差程序建立M文件夹：functiony=tongji50(A,t)%A1发病人数信息矩阵，t是年份2007,2008,2009,2010 a=size(A);B=zeros(2,12);fori=1:a(1)ifA(i,4)=tifA(i,1)=1ifA(i,5

48、)=1B(1,1)=B(1,1)+1;elseifA(i,5)=2B(1,2)=B(1,2)+1;elseifA(i,5)=3B(1,3)=B(1,3)+1;elseifA(i,5)=4B(1,4)=B(1,4)+1;elseifA(i,5)=5B(1,5)=B(1,5)+1;elseifA(i,5)=6B(1,6)=B(1,6)+1;elseifA(i,5)=7B(1,7)=B(1,7)+1;elseifA(i,5)=8B(1,8)=B(1,8)+1;elseifA(i,5)=9B(1,9)=B(1,9)+1;elseifA(i,5)=10B(1,10)=B(1,10)+1;elseifA(

49、i,5)=11B(1,11)=B(1,11)+1;elseifA(i,5)=12B(1,12)=B(1,12)+1;endelseifA(i,5)=1 B(2,1)=B(2,1)+1; elseifA(i,5)=2 B(2,2)=B(2,2)+1;elseifA(i,5)=3 B(2,3)=B(2,3)+1;elseifA(i,5)=4 B(2,4)=B(2,4)+1; elseifA(i,5)=5 B(2,5)=B(2,5)+1;elseifA(i,5)=6 B(2,6)=B(2,6)+1;elseifA(i,5)=7 B(2,7)=B(2,7)+1;elseifA(i,5)=8B(2,8)

50、=B(2,8)+1;elseifA(i,5)=9B(2,9)=B(2,9)+1;elseifA(i,5)=10 B(2,10)=B(2,10)+1;elseifA(i,5)=11B(2,11)=B(2,11)+1;elseifA(i,5)=12 B(2,12)=B(2,12)+1;endend end end y=B,%1年*12个月的男女发病人数 %发病人群信息行-6列A=xlsread('c 题数据 1.xls','sheet1','A2:F58926');%58925 a=size(A);B1=tongji50(A,2007);B2=ton

51、gji50(A,2008);B3=tongji50(A,2009);B4=tongji50(A,2010);B=B1'B2'B3'B4'B,%48*2矩阵，男女48个月的发病人数p=anova1(B),%下面做单因素方差分析附录3.统计各阶段发病人群年龄比例建立M-文件：发病人群信息A=xlsread('c题数据 1.xls','sheet1','A2:F58926');%58925行-6列a=size(A);B=zeros(1,6);fori=1:a(1) ifA(i,2)<40B(1,1)=B(1,1)+

52、1;%40岁以下病发人数elseifA(i,2)>=40&A(i,2)<50B(1,2)=B(1,2)+1;%40-50岁病发人数elseifA(i,2)>=50&A(i,2)<60B(1,3)=B(1,3)+1;%50-60岁病发人数elseifA(i,2)>=60&A(i,2)<70B(1,4)=B(1,4)+1;%60-70岁病发人数elseifA(i,2)>=70&A(i,2)<80B(1,5)=B(1,5)+1;%80岁以上病发人数elseB(1,6)=B(1,6)+1;endendB=B/a(1)% 年龄比率附录4、不同职业病发者比例%发病人群信息行