全国数学建模大赛c题.pdf

《全国数学建模大赛c题.pdf》由会员分享，可在线阅读，更多相关《全国数学建模大赛c题.pdf（31页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

2、择一项填写） ： C我们的参赛报名号为（如果赛区设置报名号的话） ：139C01所属学校（请填写完整的全名） ：浙江工贸职业技术学院参赛队员 (打印并签名) ：1.郑济明 2.王庆松_ 3.朱松祥指导教师或指导教师组负责人 ( 打印并签名 )：王积建日期： 2012年 9 月 10 日赛区评阅编号（由赛区组委会评阅前进行编号）：_20112011 高教社杯全国大学生数学建模竞赛高教社杯全国大学生数学建模竞赛编编 号号 专专 用用 页页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委

3、会评阅前进行编号）：脑卒中发病环境因素分析及干预脑卒中发病环境因素分析及干预摘摘要要_关键词关键词：_一、一、问题重述问题重述 21 世纪人类倡导人与自然和谐发展， 环境因素成为影响健康的重要因素。 脑卒中（俗称脑中风）就是与环境因素紧密相关且威胁人类生命的疾病之一。这种疾病的诱发已经被证实与环境因素有关，其中与气温和湿度存在着密切的关系。对脑卒中的发病的环境因素进行分析，其目的是为了进行疾病的风险评估， 对脑卒中高危人群能够及时采取干预措施，也让尚未得病的健康人， 或者亚健康人了解自己得脑卒中风险程度，进行自我保护。同时，通过数据模型的建立，掌握疾病发病率的规律，对于卫生行政部门和医疗机构合

4、理调配医务力量、改善就诊治疗环境、配置床位和医疗药物等都具有实际的指导意义。现从中国某城市各家医院 2007 年 1 月至 2010 年 12 月的脑卒中发病病例信息以及相应期间当地的逐日气象资料（Appendix-C2）和 数据（见Appendix-C1） 。需解决一下几个问题：问题一：根据病人基本信息，对发病人群进行统计描述。问题二：建立数学模型研究脑卒中发病率与气温、气压、相对湿度间的关系。问题二：查阅和搜集文献中有关脑卒中高危人群的重要特征和关键指标， 结合 1、2 中所得结论，对高危人群提出预警和干预的建议方案。二、问题分析脑卒中（俗称脑中风）作为威胁人类生命的疾病之一，并且病发的人

5、群受环境因素的影响不断扩展。对脑卒中人群及受环境因素的影响分析来对疾病的风险评估，对脑卒中高危人群能够及时采取干预措施成为一项无疑是一项十分复杂的系统工程。对于问题一， 利用中国某城市各家医院 2007 年 1 月至 2010 年 12 月的脑卒中发病病例信息以及相应期间当地的逐日气象资料（Appendix-C2）和数据（见Appendix-C1） 。通过 excel 对已知数据进行统计整理，再利用 matlab 程序对脑卒中病发者的性别指数、年龄指数、职业指数、月份指数进行合理的统计得出相应数据比率。三、模型假设三、模型假设3.13.1 模型假设：1) 发病病例的信息中，若两个病例的信息相同

6、，则视为不同的两个人；2）以诊断报告时间为准来统计发病人群的数量；_3）导致脑卒中发病的内在原因只与性别、年龄、职业有关；4）导致脑卒中发病的外在原因只与气压、温度和湿度有关；5）气压、温度和湿度之间具有相关关系；6）月平均气压、月平均最高气压、月平均最低气压具有相关关系；7）月平均温度、月平均最高温度、月平均最低温度具有相关关系；8）月平均湿度、月平均最高湿度、月平均最低湿度具有相关关系；9）关于环境因素如气压、温度和湿度的观测数据都是准确可靠的；10）按照国际惯例1,发病率以 10 万人群的发病人数来表示。但由于本题是研究某地区的发病人数，并没有与其它地区比较，所以在本题分析中，发病率以发

7、病人数来表示。四、符号说明四、符号说明定义 1，月平均气压是日平均气压的平均值。月平均最高气压是日平均最高气压的平均值。月平均最高气压是日平均最高气压的平均值。定义 2，月平均温度是日平均温度的平均值。月平均最高温度是日平均最高温度的平均值。月平均最高温度是日平均最高温度的平均值。定义 3，月平均湿度是日平均湿度的平均值。月平均最低湿度是日平均最低气压的平均值。N1表示男性病例总数，N2表示女性病例总数，N表示总病例数;五、模型的建立及求解五、模型的建立及求解5.05.0 发病人群数据的预处理发病人群数据的预处理根据已知题意给出的中国某城市各家医院 2007 年 1 月至 2010 年 12

8、月的脑卒中发病病例信息以及相应期间当地逐日气象资料，进行如下数据预处理： 1）以 2007 年 1 月至 2010 年 12 月的一共 48 个月的脑卒中发病病例为准，其他时间数据应当删除，一共得到 58925 个病例.2）如果病例的信息中，年龄与职业不符（例如：12 周岁是老师） 、诊断时间不详、数据明显出错的都不应该考虑在统计范围之内，应当删掉。3）5.15.1 对发病人群的统计分析对发病人群的统计分析（问题 1）5.1.15.1.1 性别分析1 1）性别差异性简单分析）性别差异性简单分析男、女性病发比例为Nx1i1i,i 1,2（1）N其 中 ，i 1表 示 男 性 ，i 2表 示 女

9、性 。 经 统 计 ，N 58925，N11 31832，N12 27093，代入（1）得男、女病发比例分别为 54.02%和 45.98%（matlab 程序见附录 1） 。可见男性在脑卒中的病发者要大于女性脑卒中病_发人数。2)2)单因素方差分析单因素方差分析22逐月统计男女病例人数，考察在相同时间点上男女人群发病人数是否有显著差异，给定显著性水平 0.05，分析结果为F 5.54，对应的p 0.0206 0.05（见图1 ）,又 查 表 得F2(r 1,n r) F0.05(2 1,48 2) F0.05(1,46) 4.08, 由 于F F0.05(1,46),所以脑卒中发病男女人群有

10、显著差异(matlab 程序见附录 2)。图 1 男女发病人群的单因素方差分析结果5.1.2 不同年龄段发病人群差异性分析1 1）简单分析）简单分析不同年龄阶段发病比例为Nx2i2i,i 1,2,3,4,5,6（2）Ni 1,2,3,4,5,6分别表示其中，“40 岁以下” 、“40-50” 、“50-60” 、“60-70” 、“70-80” 、“80 以上” 。经统计，将不同年龄阶段脑卒中病发者人数代入（ 2）式，得到不同年龄阶段脑卒中病发者比例，见图 2。 （matlab 程序见附录 3）40%30%20%10%0%40以下40-5050-6060-7070-8080以上34.06%23

11、.19%12.88%1.75%4.73%23.39%图 2 不同年龄段发病人群比例图由图 2 可以看出在 50 岁以下的人口中脑卒中病发的人数比例较小， 70-80 之间脑卒中的比例最为严重， 80 岁以上的人脑卒中较为严重， 所以高龄的人是发生脑卒中的高危人群，我们应当高度关注。2)2)单因素方差分析单因素方差分析根据图 2 结果，剔除“40 岁以下”和“40-50”年龄段，对其余 4 个年龄段进行单因素方差分析，逐月统计不同年龄段发病人群人数， 考察在相同时间点上不同年龄_段发病人群人数是否有显著差异，给定显著性水平 0.05，分析结果为F 45.6，对应的p 0.0000 0.05（见图

12、 3） ,所以脑卒中不同年龄段发病人群有显著差异(matlab 程序见附录 4（tongji6.m 和 tongji60.m）)。图 3 不同年龄段发病人群的单因素方差分析结果 5.1.3 不同职业发病人群的差异性分析1 1）简单分析）简单分析不同职业发病比例为Nx3i3i,i 1,2,.,9（3）N其中，i 1,2,3,4,5,6,7,8,9分别表示“农民” 、 “工人” 、 “退休人员” 、 “教师” 、 “渔民” 、 “医务人员” 、 “职工” 、 “离退人员” 、 “其它职业” 。经统计，不同职业脑卒中病发者的比例，见图 4。 （matlab 程序见附录 4）60%48.06%50%4

13、0%29.37%30%20%7.28%10.70%10%0.36%0.10%0.14%1.19%2.80%0%农民工人教师渔民退休人员医务人员职工离退人员其他职业图 4不同职业病发者比例由图 4 得出农民、工人、退休人员、其他职业的人员患脑卒中的比例偏高，说明了职业也是患脑卒中的重要因素。2)2)单因素方差分析单因素方差分析根据图 4 结果，对农民、工人、退休人员进行单因素方差分析，给定显著性水平 0.05，分析结果为F 95.36，对应的p 0.0000 0.05（见图 5） ,所以脑卒中不同年龄段发病人群有显著差异(matlab 程序见附录 5（tongji7.m 和 tongji70.m

14、）)。_图 5 不同年龄段发病人群的单因素方差分析结果5.1.4 不同月份发病者的差异性分析1 1）简单分析）简单分析定义季节指数为Sx4i4i,i 1,2,.,12（4）S其中，S4i为第i月的平均人数，S为 48 个月的月平均人数。经统计，不同月份脑卒中病发者的比例，见图 6。 （matlab 程序见附录 6）1.210.80.60.40.201月2月3月4月5月6月7月8月9月10月11月12月系列10.70930.88041.01440.99181.10121.03381.11891.03461.03151.03331.00991.0409图 62007 年-2010 年各月季节指数由

15、图 6 看出在五、六、七月份为脑卒中高发期，一、二月为低发期。5.2 发病率与气压、气温、相对湿度间的关系分析（问题 2）由于题目提供了环境因素（气压、温度和湿度）的8 个变量，根据假设5）8） ，这 8 个变量间具有明显的显著相关关系，所以必须做降维处理，把8 个变量整合成互不相关的少数几个变量，然后再寻找发病率与这少数几个变量的关系式。 这需要进行主成分分析。5.2.1 主成分分析法31 1）基本原理基本原理主成分分析是把多个变量转化为少数几个新综合变量的一种多元统计方法，其基本思想就是在保留原始变量尽可能多的信息的前提下达到降维的目的， 从而简化问题的复杂性并抓住问题的主要矛盾.其手段是

16、将原来众多的具有一定相关性的变量重新组合成新的少数几个相互无关的综合变量（也叫抽象变量） ，来代替原来变量，这些新的综合变量称之为主成分.一般地说， 利用主成分分析得到的主成分与原来的变量之间有如下基本关系： （1）每一个主成分都是各原始变量的线性组合.（2）主成分的数目大大少于原始变量的数目.（3）主成分保留了原始变量的绝大多数信息.（4）主成分之间互不相关.据此我们建立数学模型.2)2)数学模型数学模型在一个统计问题中，假设我们收集到n个样品，每个样品观测到p个变量（记为x1,x2,xp,为简单起见，可以设xi均值为 0，方差为1，(1i p） ，构成一个n p阶的样本原始资料阵X xij

17、np._主成分分析的目的在于利用p个原始变量（x1,x2,xp）构造少数几个新的综合变量，使得新变量为原始变量的线性组合，新变量互不相关，新变量包含p个原始变量的绝大部分信息.这样定义x1,x2,xp为原始变量，y1, y2, ym(m p)为新的综合变量指标，每一个新综合变量指标是p个原始变量的线性组合：y1 a11x1a12x2a1pxpy2 a21x1a22x2a2pxp (5)y a x ax axm1 1m22mppm同时要求满足以下几个条件： （1）yi与yj相互无关； （2）y1是x1,x2,线性组合中方差最大者；y2是y1与不相关的x1,x2, ym是z1,z2,zm1分别都不

18、相关的x1,x2,xp的一切,xp的所有线性组合中方差最大者；,xp的所有线性组合中方差最大者.则新变量y1, y2, ym分别称为原变量x1,x2,xp的第一、第二、第m主成分., p)在从以上的分析可以看出，主成分分析的实质就是确定原来变量xj( j 1,2,诸主成分yi(i 1,2,m)上的系数aij(i 1,2,m; j 1,2, p).从数学上可以证明，他们分别是p个原始变量（x1,x2,xp）相关矩阵的前m个具有较大特征值所对应的特征向量，而各个新综合变量yi的方差var(yi)恰好是相应的特征值i.各主成分的方差贡献大小按特征根顺序排列， 是依次递减的， 即12p 0.其几何意义

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

20、主成分分析.Bartlett 球形检验是检验相关矩阵是否是单位矩阵，即各变量是否各自独立._（3）对原始数据进行标准化，即令*xijxij xjsj（6）其中xj,sj分别为第j列元素的样本均值和样本标准差，即1n1nxjxij,sj(xij xj)2ni1n1i1*则X* (xij)np为标准化的样本资料库.（4）由标准化后的数据矩阵求协方差矩阵，或者由原始数据矩阵求相关系数矩阵R.这两种方法结果相等.本文采用直接计算原始数据的相关矩阵的方法 （对于数量级差别较大或者有量纲的数据宜适用）.设原始数据X的相关系数矩阵为r11r12.r1prr.r21222p（7）R .rr.rnpn1n2ri

21、j(i, j 1,2, p)为原变量xi与xj的相关系数，rij rji，其计算公式为rij(xk1nki xi)(xkj xj)n(xk1n（8）ki xi)2(xkj xj)2k1（5）计算R的特征根和特征向量；根据特征方程E R 0得R的特征根为i(i 1,2,., p)，将特征根按照从大到小的顺序排列，排列后的特征根不妨仍然表示为12.p 0.同时可得对应的特征向量u1,u2,.,up，将他们标准正交化，u1,u2,.,up称为主轴.（6）计算所有变量的方差贡献率及累计方差贡献率；i的方差贡献率为eiii1pi 1,2,i, p（9）i的累计方差贡献率为_Eik1pi1iki 1,2,

22、 p, k 1,2,.,m,m p（10）i（7） 确定主成分的数目m. 方法有： 一般取累计贡献率达 85%95%的主成分；选用所有i1的主成分； 累计特征值乘积大于 1 的主成分； 画出特征值变化曲线，以转折点位置为标准判断.本文采用累计贡献率达 85%95%的主成分.（8）确定主成分函数表达式模型. 设m个主成分对应的特征向量分别为A1、A2、.Am，其中Aj a1ja2 j.ap j，ak j表示Aj的第k行的元素，则第j个T主成分yj的函数表达式为 x1x2TyjAj a . 1jxpa2 j x1px2.ap jak jxk（11） . k1xp（9）提炼主成分yj的抽象意义.由x

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

24、的平均值、 月最高气温的平均值、 月平均气压的平均值 8 项指标，并分别记为x1，x2x8. 每个指标有 48 个数据（见附件 1） 。使用 SPSS 软件进行求解（见附录 7） 。（2）将原始数据标准化， （SPSS 内部计算）.（3）求原始数据的相关系数矩阵R，如图 7 所示._图 7 相关系数矩阵图 8 因子分析检验图从图 8 看出，表格的第一行为检验变量间偏相关程度的 KMO 统计量，其值在 0.6之上才适合做主成分分析，效果显著，如果小于 0.6，效果不显著，不适合做主成分分析。下面的三行为球形检验的结果，球形检验原假设的变量是不相关的，显然只有拒绝原假设的情况下数据才适合做因子分析

25、。 本例中 KMO 值为 0.720， 球形检验显著，两个条件都满足，变量间相关程度大，适合做因子分析。（4）计算矩阵R的特征根、各因子的方差贡献率及累计方差贡献率，并确定主成分的个数.如图 9 所示。图 9R特征值及其累计方差贡献率从图 9 中可以看出，第一、第二主成分对方差的累计贡献率达到 95.461%，它们分别对应着原样本数据点数据变异的最大、 次大方向， 是原变量系统的一个最佳整合，从而我们可以以 95.461%的精度将变量的有效维数从 8 维降至 2 维.因此可以将前 2个因子作为主因子.（5）确定主成分函数表达式模型，因子得分系数矩阵如图 10 所示._图 10因子得分系数矩阵设

26、 2 个主成分分别为y1, y2，则建立模型为y1 0.184x10.183x20.185x30.175x40.186x50.164x60.096x70.121x8（12）y2 0.048x10.042x20.057x30.020 x40.055x50.015x60.536x70.567x8其中x1,x2,x8.均为原变量经过均值为 0，方差为 1 标准化后的变量.（6）对主成分y1, y2的意义进行解释。图 11 给出了原变量与第 1、第2 主成分的相关系数.图 11 旋转后的因子载荷矩阵第一主成分y1， 与原变量x1(平均气压的平均值)、x2（最高气压的平均值）x3（最低气压的平均值） 、

27、x6（最低气温的平均值）的相关系数的绝对值都超过了 0.948，因此它是一个反映气温和气压的综合因子，我们称之为气压温度因子.第二主成分y2，与原变量x7（月平均相对湿度的平均值）的相关系数为 0.925、x8（月最低相对湿度的平均值）的相关系数为 0.948，其余的都不超过 0.266，因此它是一个反映相对湿度的因子，称为湿度因子.（7）计算 2 个主成分的函数值.将 48 个经过标准化的数据xi*代入模型yj， 可以得到 48 个地区的主成分yj的函数_值，结果如表 1 所示.表 1 48 个月对应的主成分函数值序号1234567891011121314151617181920212223

28、24y1-1.6817y20.5485z48752878485810821137104310801206126512021297120917161787159817551533207813221384133114751472序号252627282930313233343536373839404142434445464748y1-1.121y21.77z94679973986185578875296282876966762384112801671155317151618162117161647170916151719-1.601-0.061-0.6750.00760.1144-0.9170.92

29、19-2.0981.34170.05841.35210.47150.94021.66440.35751.51150.2595-0.727-1.0911.1667-1.3310.3153-1.4530.2577-1.0830.8107-0.8230.1507-0.401-0.0760.7358-0.6590.78240.70991.24550.64451.3537-0.1330.58160.9467-0.2460.4312-0.531-0.633-0.593-1.472-0.7683-0.2082-0.4507-0.21660.1520.92591.52390.5675-0.1318-1.189

30、3-1.59090.75540.20971.00550.1920.53041.1598-1.77581.5284-0.4272-0.6524-0.9464-1.5177-0.4902-1.3599-1.3556-0.1654-2.01870.18450.67641.14590.5359-0.1899-0.9077-1.031.86590.73561.14160.91490.57660.9955-1.31491.6182-0.2525-1.036-0.9904（8）检验主成分模型.由于主成分分析的 4 个条件中的前 3 个（每一个主成分都是各原始变量的线性组合；主成分的数目大大少于原始变量的数目

31、；主成分保留了原始变量的绝大多数信息） ，只要检验 4 个主成分是否相关即可.由步骤（7）计算的 2 个主成分的得分矩阵Y yij482，求矩阵Y的协方差矩阵如图 12 所示._图 12 因子得分的协方差矩阵从图 12 可以看出，主成分得分的协方差矩阵为单位矩阵，说明提取的 2 个主成分是互不相关的.满足假设的条件，模型和结果有效.5.2.2 多元非线性回归分析4将 48 个月的发病率作为因变量，记作z，将发病率的 48 个数据填入表 1 中。下面寻找发病率z与主成分y1, y2的关系式，这需要使用多元非线性回归分析方法。经过反复试验探索，找到的非线性回归模型为a ya4y2az expa1y

32、1a2y2315a6（13）y1 y2y1(1 y2)y1其中，a1 0.1281,a2 0.0472,a3 0.0273，a4 0.0021,a5 0.0407,a6 7.0361。模型检验的p 0.0399 0.05，说明模型有效。平均绝对相对误差为 26.38%。5.2.3 结果分析1）从非线性回归模型（12）可以得到以下结论：（1）由a1 0可得，发病率与气压温度因子具有正相关性；（2）由a2 0可得，发病率与湿度因子具有负相关性；（3）由a1 a2可得，气压温度因子比湿度因子对于发病率的影响显著；2）从主成分模型（11）可以得到以下结论：由第 1 个方程可知：（1）由x1,x2,x3

33、的系数为负值可得，气压温度因子与月平均气压、月平均最高气压、月平均最低气压具有负相关性；（2）由x4,x5,x6的系数为正值可得，气压温度因子与月平均温度、月平均最高温度、月平均最低温度具有正相关性；（3）由x7,x8的系数为负值可得，气压温度因子与月平均相对湿度、月平均最低相对湿度具有负相关性；由第 2 个方程可知：（4）由x1,x2,x3的系数为正值可得，湿度因子与月平均气压、月平均最高气压、月平均最低气压具有正相关性；（5）由x4 0,x5 0,x6 0可得，湿度因子与月平均温度、月平均最高温度成负相关性，与月平均最低温度具有正相关性；（6）由x7,x8的系数为正值可得，湿度因子与月平均

34、相对湿度、月平均最低相对湿度具有正相关性；_5.3 高危人群预警分析（问题 3）首先预测 2011 年的气象状况， 然后预测未来 2011 年的高危人群的发病率和发病时间，最后提出预警和干预措施。5.3.1 未来 2011 年气象预测1 1）气压、温度状况预测）气压、温度状况预测观察月平均气压x1、月平均最高气压x2、月平均最低气压x3、月平均温度x4、月平均最高温度x5，月平均最低温度x6的历史数据，发现随时间做周期性变化，于是建立余弦函数模型x Acos(t ) B， 利用过去 36 个月的历史数据进行参数估计， 然后使用2009年的12个数据进行预测， 评估误差并检验模型的可靠性。 最后

35、预测出2010年的 12 个月的数据。建模结果见表 2.（MATLAB 程序见附录 8） 。2011 年 12 个月的预测结果见表 3.表 2 模型参数估计结果变量表达式参数拟合值平均相对误差0.15%x1x Acos(t ) BA 10.6, 5, 1,B 1015.5A 11,2, 5, 4,B 1017.8A 10,1, 5, 3,B 1013.2A 12.1375, 0.5211, 2.1851,B 16.8618x2x Acos(t ) B0.14%x3x Acos(t ) B0.15%x4x Acos(t ) B11.77%x5A 12.0959, 0.5198,x Acos(t

36、) B 6.8425,B 21.0561A 12.2819, 0.5224,x Acos(t ) B 0.6683,B 13.54827.75%x634.71%_表 3 2011 年 12 个月的气压、温度预测值月份平均温度平均压力最高压力最低压力最高温度最低温度月份平均温度平均压力最高压力最低压力最高温度最低温度14.99781026.11018.21014.723.742616.0265728.68621004.91017.61011.718.098510.983925.29691024.31023.71019.517.53029.6937828.48031006.71012.11006.

37、924.314417.319138.66621020.11027.71022.612.2494.3892925.19011010.91008.11003.829.669722.6484414.21111014.710291023.29.2941.52791019.6891016.21006.61003.232.7525.5502520.45961009.61027.41021.19.44561.87311113.43741021.41008.21005.232.741725.2503625.753110061023.21016.912.66385.3328128.09481025.11012.

38、31009.429.64721.82892 2）湿度状况预测）湿度状况预测观察月平均湿度x7、 月平均最低湿度的历史数据x8， 发现它们随时间做平稳性波动，于是建立马尔克夫模型，利用过去 36 个月的历史数据进行建模，然后预测出第37 个数据。采用“新陈代谢”思想，把第 37 个数据加入建模序列，并同时去掉最老的第 1 个数据，保持数据“等维” ，建模并预测出第 38 个数据，如此滚动预测，直至预测出 2009 年的 12 个数据，并做误差分析，检验模型的可靠性。最后预测出 2011年的 12 个数据。（1）自相关系数原始序列X(0)=X(0)(1),X(0)(2),., X(0)(n)的各阶

39、自相关系数反映已知数据对未来数据的影响程度. 各阶自相关系数为nwrwXk1(0)(0)(0)(k) X(0) X(k w) XXk1n(0)(k) X(0)2（13）式中，X(0)1n(0)=X(k)（14）nk1对各阶自相关系数归一化得，w=rwrw1t,w 1,2,.,t（15）w_w可作为各阶步长的马尔柯夫链权重，t是按预测需要计算的最大阶数，一般取rw 0.3.根据rw 0.3可以确定转移步数w.（2 2）加权马尔柯夫模型）加权马尔柯夫模型状态划分。设划分的m个湿度区间为E1a11,a21,Eia1i,a2i,i 2,3,.,m 1Ema1m,a2m其中，a11尽可能小，a2m尽可能

40、大.，如果(k)Ei,i 1,2,., m则表明第k年的相对误差处于第i种状态.(w)状态转移概率矩阵的构造。设w步转移概率为pij，记：(w)mij(w)pij=Mi,i, j 1,2,.,m（15）(w)其中，mij表示状态Ei经过w步转移到状态Ej的次数，Mi为状态Ei出现的次数.由于数据序列最后的状态转向不确定，故计数Mi时要去掉数据序列中最末的w个数据(也就是只考虑前面的n w个数据).(w)由pij构成的矩阵称为w步转移概率矩阵，记作w)p1(m(w)p2m（16）R(w)(w)(w)pm.p2mm已知每一步的概率转移矩阵和每一步的初始状态，则马尔柯夫链就可以确定.预测值计算(w)

41、 p11(w)p21.(w)pm1(w)p12.(w)p22.选取距离预测年最近的t(t m)个年份， 按照距离预测年由近到远， 转移步数w分别为1,2,.,t，以这几年的相对误差所对应的状态为初始状态，不妨设第1,2,.,t年所(1)(2),E2,., Et(t)，其中，w1,2,.,m.例如，当2 5时，对应的初始状态分别为E1(2)(2)E2=E5，说明距离预测年第 2 年的状态是第 5 状态.在转移步数w对应的转移矩阵(w)(w)(w)(w)(w)R(w)中，取起始状态Ew.，pwm，从而组成新的概所对应的行向量pw=pw1，pw2，_率矩阵(1) p11(2)pR 21.(t)pt1

42、(1)p12.(2)p22.t)pt(2.(1)p1m(2)p2 m（17）)pt(tm将矩阵R加权得(1)(1)(1)1p111p12.1p1m(2)(2)(2)pp.p22222 m（18）R221.(t)(t)(t)pp.ptt2ttmtt1将矩阵R按列求和得ttt(w)(w)(w)p p1, p2,., pm=wpw1,wpw2,.,wpwm（19）w1w1w1找出向量p的最大分量得pM maxp1, p2,., pm,M 1,2,.,m（20）分量pM所对应的状态EM就是预测年的状态，则该年度的预测值为Ma1Ma2M（21）2（3）计算过程和结果以预测 2010 年第 1 月的数据为

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

44、步和 2 步内的转移概率矩阵分别见表 5 和m表 5 1 步转移概率矩阵R(1)2/301/80001/31/51/41/70002/51/83/71/91/301/501/74/91/301/53/82/71/30001/801/91/3表 6 2 步转移概率矩阵R(2)1/31/51/80001/301/82/71/801/3004/71/21/301/55/800002/51/82/71/42/301/501/71/80组成预测年份的新转移概率矩阵.选择离预测年最近的 2 个年份，转移步数分别为w1,2，根据（17）式得预测年的转移概率矩阵R37，见表 7.表 7月份状态步长363523

45、12权重10.66250.3375000预测年的转移概率矩阵R3721/71/833/71/242/7051/71/4601/8概率来源R(1)R(2)加权求和0.13680.45270.09460.27370.0422确定预测年份的状态确定预测年份的状态. .预测年的状态向量的最大分量值为 0.4527，对应的状态为第 3 状态，即第 37 个月的湿度将处于第 3 状态，湿度 67.5.将第 37 个月的湿度值 67.5 放入序列中，同时去掉第1 个月的湿度数据，重新构建马尔柯夫链，得第 38 个月的湿度。以此类推，可得 2010 年 12 个月的湿度值，见表 8.平均绝对相对误差为 7.9

46、7%，可靠性高.表 8 2010 年湿度预测月份实际x7167.5277.5372.5467.5577.5672.5767.5877.5972.51067.5111277.5725_70.944.85%74.683.78%70.812.39%69.572.98%68.3913.32%75.5810.69%71.588.27%74.299.14%62.2316.50%7模拟x绝对相对误差76.45.10%77.26.09%68.912.48%表 8 2011 年湿度预测月份123456789101112x772.577.572.567.577.572.567.577.572.567.577.57

47、2.5x857.557.557.557.557.557.557.557.557.557.557.557.55.3.2 未来 2011 年 12 个月的发病人数预测将预测得到的 2011 年 12 个月的气象数据代入主成分模型（12）和多元非线性回归模型（13） ，计算得到 2011 年发病率的发展趋势，如图 13 所示。0.60.50.40.30.20.101234567891011122011年月份发病率图 132011 年发病率趋势从图 13 可知，2011 年发病率的发展趋势是，15 月保持平稳态势，从6 月份开始呈现上升态势，至 8、9、10 月份到达高发期，之后，11 和 12 月份

48、下降进入低发期。5.3.3 干预措施和建议查阅文献5-9，发现以下信息：1）脑卒中,又称中风或脑血管意外,它包括脑出血、蛛网膜下腔出血、脑梗死和短暂性脑缺血发作等急性脑血管病 ,是一组突然起病,以出现意识障碍和局灶性神经功能缺失为共同特征的急性脑血管病。2）它具有以下特征:(1)发病率高;(2)致残率高;(3)死亡率高;(4)复发率高。3）它的关键指标为： （1）头晕。 （2）肢体麻木。 （3）暂时性吐字不清或讲话不灵。 （4）肢体无力或活动不灵。 （5）与平时不同的头痛。 （6）不明原因突然跌倒或晕倒。 （7）短暂意识丧失或个性和智力的突然变化。 （8）全身明显乏力，肢体软弱无力。（9）恶心

49、呕吐或血压波动。脑卒中病人大部分是“三高” 、有家族病史与前科老年人。建议措施： （1）卫生机构在 2011 年的秋季前做好迎接高发病人的准备。 （2）预防为主，多锻炼，多吃蔬菜，不吸烟，不喝酒，防“三高” 。六、需要进一步研究的问题六、需要进一步研究的问题对高危人群如男性人群、60 岁以上人群、农民群体进行预测，做好预警和干预。_七、模型评价七、模型评价7.17.1 模型优点：1）主要因素并进行了定量分析，7.2 模型缺点：1)当考虑定性相关因素较多时，八、参考文献1周晓平，杨进.脑卒中发生时气节规律及其气象医学原理探讨.中医杂志，2006 年 6 月,第 47 卷第 6 期.2谢文龙,尚涛

50、.SPSS 统计分析与数据挖掘.北京：北京电子工业出版社，2012.13韩中庚.数学建模方法与应用.北京:高等教育出版社，2005.6(2006 重印).4陈在余,陶应虎.统计学原理与实务.北京:清华大学出版社，2009.3复发高危人群，低温危害高温危害九、附录附录 1.统计发病人群的性别比例建立 M-文件：%统计发病人群的性别比例A=xlsread(c 题数据 1.xls,sheet1,A2:F58926);%读入发病人群信息，58925行-6 列a=size(A);B=zeros(1,2);for i=1:a(1) if A(i,1)=1 B(1,1)=B(1,1)+1;%男性病人的人数