全国数学建模大赛2012C题

2011高教社杯全国大学生数学建模竞赛

承 诺 书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C 我们的参赛报名号为（如果赛区设置报名号的话）： 139C01 所属学校（请填写完整的全名）： 浙江工贸职业技术学院 参赛队员 (打印并签名) ：1. 郑济明 2. 王庆松 3. 朱松祥 指导教师或指导教师组负责人 (打印并签名)： 王积建

日期： 2012 年 9 月 10 日

赛区评阅编号（由赛区组委会评阅前进行编号）：

2011高教社杯全国大学生数学建模竞赛

编 号 专 用 页

评 阅 人 评 分 备 注 赛区评阅编号（由赛区组委会评阅前进行编号）：

赛区评阅记录（可供赛区评阅时使用）：

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

脑卒中发病环境因素分析及干预 摘 要

关键词：

一、问题重述

21世纪人类倡导人与自然和谐发展，环境因素成为影响健康的重要因素。脑卒中（俗称脑中风）就是与环境因素紧密相关且威胁人类生命的疾病之一。这种疾病的诱发已经被证实与环境因素有关，其中与气温和湿度存在着密切的关系。对脑卒中的发病的环境因素进行分析，其目的是为了进行疾病的风险评估，对脑卒中高危人群能够及时采取干预措施，也让尚未得病的健康人，或者亚健康人了解自己得脑卒中风险程度，进行自我保护。同时，通过数据模型的建立，掌握疾病发病率的规律，对于卫生行政部门和医疗机构合理调配医务力量、改善就诊治疗环境、配置床位和医疗药物等都具有实际的指导意义。

现从中国某城市各家医院2007年1月至2010年12月的脑卒中发病病例信息以及相应期间当地的逐日气象资料（Appendix-C2）和 数据（见Appendix-C1）。需解决一下几个问题：

问题一：根据病人基本信息，对发病人群进行统计描述。

问题二：建立数学模型研究脑卒中发病率与气温、气压、相对湿度间的关系。 问题二 ：查阅和搜集文献中有关脑卒中高危人群的重要特征和关键指标，结合1、2中所得结论，对高危人群提出预警和干预的建议方案。

二、问题分析

脑卒中（俗称脑中风）作为威胁人类生命的疾病之一，并且病发的人群受环境因素的影响不断扩展。对脑卒中人群及受环境因素的影响分析来对疾病的风险评估，对脑卒中高危人群能够及时采取干预措施成为一项无疑是一项十分复杂的系统工程。

对于问题一，利用中国某城市各家医院2007年1月至2010年12月的脑卒中发病病例信息以及相应期间当地的逐日气象资料（Appendix-C2）和数据（见Appendix-C1）。通过excel对已知数据进行统计整理，再利用matlab程序对脑卒中病发者的性别指数、年龄指数、职业指数、月份指数进行合理的统计得出相应数据比率。

三、模型假设

3.1模型假设：

1) 发病病例的信息中，若两个病例的信息相同，则视为不同的两个人； 2）以诊断报告时间为准来统计发病人群的数量；

4

3）导致脑卒中发病的内在原因只与性别、年龄、职业有关； 4）导致脑卒中发病的外在原因只与气压、温度和湿度有关； 5）气压、温度和湿度之间具有相关关系；

6）月平均气压、月平均最高气压、月平均最低气压具有相关关系； 7）月平均温度、月平均最高温度、月平均最低温度具有相关关系； 8）月平均湿度、月平均最高湿度、月平均最低湿度具有相关关系； 9）关于环境因素如气压、温度和湿度的观测数据都是准确可靠的；

10）按照国际惯例[1],发病率以10万人群的发病人数来表示。但由于本题是研究某地区的发病人数，并没有与其它地区比较，所以在本题分析中，发病率以发病人数来表示。

四、符号说明

定义1，月平均气压是日平均气压的平均值。月平均最高气压是日平均最高气压的平均值。月平均最高气压是日平均最高气压的平均值。

定义2，月平均温度是日平均温度的平均值。月平均最高温度是日平均最高温度的平均值。月平均最高温度是日平均最高温度的平均值。

定义3，月平均湿度是日平均湿度的平均值。月平均最低湿度是日平均最低气压的平均值。

N1表示男性病例总数，N2表示女性病例总数，N表示总病例数;

五、模型的建立及求解

5.0发病人群数据的预处理

根据已知题意给出的中国某城市各家医院2007年1月至2010年12月的脑卒中发病病例信息以及相应期间当地逐日气象资料，进行如下数据预处理：

1）以2007年1月至2010年12月的一共48个月的脑卒中发病病例为准，其他时间数据应当删除，一共得到58925个病例.

2）如果病例的信息中，年龄与职业不符（例如：12周岁是老师）、诊断时间不详、数据明显出错的都不应该考虑在统计范围之内，应当删掉。 3）

5.1对发病人群的统计分析（问题1） 5.1.1性别分析

1）性别差异性简单分析

男、女性病发比例为

Nx1i?1i,i?1,2 （1）

N其中，i?1表示男性，，i?2表示女性。经统计，

N?58925，N11?31832，N12? 27093，代入（1）得男、女病发比例分别为54.02%和45.98%（matlab程序见附录1）。可见男性在脑卒中的病发者要大于女性脑卒中病

5

发人数。

2)单因素方差分析[2]

逐月统计男女病例人数，考察在相同时间点上男女人群发病人数是否有显著差异，给定显著性水平??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),所以脑卒中发病男女人群有显著差异(matlab程序见附录2)。

图1 男女发病人群的单因素方差分析结果

5.1.2不同年龄段发病人群差异性分析

1）简单分析

不同年龄阶段发病比例为

Nx2i?2i,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）

400 %0@以下40-5050-6060-7070-8080以上

34.06#.19.88%1.75%4.73#.39%

图2 不同年龄段发病人群比例图

由图2可以看出在50岁以下的人口中脑卒中病发的人数比例较小，70-80之间脑卒中的比例最为严重，80岁以上的人脑卒中较为严重，所以高龄的人是发生脑卒中的高危人群，我们应当高度关注。

2)单因素方差分析

根据图2结果，剔除“40岁以下”和“40-50”年龄段，对其余4个年龄段进行单因素方差分析，逐月统计不同年龄段发病人群人数，考察在相同时间点上不同年龄

6

段发病人群人数是否有显著差异，给定显著性水平??0.05，分析结果为F?45.6，对应的p?0.0000?0.05（见图3） ,所以脑卒中不同年龄段发病人群有显著差异(matlab程序见附录4（tongji6.m和tongji60.m）)。

图3 不同年龄段发病人群的单因素方差分析结果 5.1.3不同职业发病人群的差异性分析

1）简单分析

不同职业发病比例为

Nx3i?3i,i?1,2,...,9 （3）

N其中，i?1,2,3,4,5,6,7,8,9分别表示“农民”、“工人”、“退休人员”、“教师”、“渔民”、“医务人员”、“职工”、“离退人员”、“其它职业”。经统计，不同职业脑卒中病发者的比例，见图4。（matlab程序见附录4）

60H.06P@).370 %7.28.70%0.36%0.10%0.14%1.19%2.80%0%农民工人教师渔民退休人员医务人员职工离退人员其他职业

图4 不同职业病发者比例

由图4得出农民、工人、退休人员、其他职业的人员患脑卒中的比例偏高，说明了职业也是患脑卒中的重要因素。

2)单因素方差分析

根据图4结果，对农民、工人、退休人员进行单因素方差分析，给定显著性水平

??0.05，分析结果为F?95.36，对应的p?0.0000?0.05（见图5） ,所以脑卒中不同年龄段发病人群有显著差异(matlab程序见附录5（tongji7.m和tongji70.m）)。

7

图5 不同年龄段发病人群的单因素方差分析结果

5.1.4不同月份发病者的差异性分析

1）简单分析 定义季节指数为

Sx4i?4i,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 图6 2007年-2010年各月季节指数

由图6看出在五、六、七月份为脑卒中高发期，一、二月为低发期。 5.2发病率与气压、气温、相对湿度间的关系分析（问题2）

由于题目提供了环境因素（气压、温度和湿度）的8个变量，根据假设5）～8），这8个变量间具有明显的显著相关关系，所以必须做降维处理，把8个变量整合成互不相关的少数几个变量，然后再寻找发病率与这少数几个变量的关系式。这需要进行主成分分析。

5.2.1主成分分析法[3] 1）基本原理

主成分分析是把多个变量转化为少数几个新综合变量的一种多元统计方法，其基本思想就是在保留原始变量尽可能多的信息的前提下达到降维的目的，从而简化问题的复杂性并抓住问题的主要矛盾.其手段是将原来众多的具有一定相关性的变量重新组合成新的少数几个相互无关的综合变量（也叫抽象变量），来代替原来变量，这些新的综合变量称之为主成分.

一般地说，利用主成分分析得到的主成分与原来的变量之间有如下基本关系：（1）每一个主成分都是各原始变量的线性组合.（2）主成分的数目大大少于原始变量的数目.（3）主成分保留了原始变量的绝大多数信息.（4）主成分之间互不相关.据此我们建立数学模型. 2)数学模型

在一个统计问题中，假设我们收集到n个样品，每个样品观测到p个变量（记为，构成一个n?p阶x1,x2,?xp,为简单起见，可以设xi均值为0，方差为1，(1?i?p）的样本原始资料阵X??xij?n?p.

8

主成分分析的目的在于利用p个原始变量（x1,x2,?,xp）构造少数几个新的综合变量，使得新变量为原始变量的线性组合，新变量互不相关，新变量包含p个原始变量的绝大部分信息.这样定义x1,x2,?,xp为原始变量，y1,y2,?,ym(m?p)为新的综合变量指标，每一个新综合变量指标是p个原始变量的线性组合：

?y1?a11x1?a12x2???a1pxp??y2?a21x1?a22x2???a2pxp ? (5)

???y?ax?ax???axm11m22mpp?m同时要求满足以下几个条件：（1）yi与yj相互无关；（2）y1是x1,x2,?,xp的一切线性组合中方差最大者；y2是y1与不相关的x1,x2,?,xp的所有线性组合中方差最大者；

?,ym是z1,z2,?,zm?1分别都不相关的x1,x2,?,xp的所有线性组合中方差最大者.则新变量y1,y2,?,ym分别称为原变量x1,x2,?,xp的第一、第二、?、第m主成分.

从以上的分析可以看出，主成分分析的实质就是确定原来变量xj(j?1,2,?,p)在诸主成分yi(i?1,2,?,m)上的系数aij(i?1,2,?,m;j?1,2,?,p).从数学上可以证明，他们分别是p个原始变量（x1,x2,?,xp）相关矩阵的前m个具有较大特征值所对应的特征向量，而各个新综合变量yi的方差var(yi)恰好是相应的特征值?i.各主成分的方差贡献大小按特征根顺序排列，是依次递减的，即?1??2????p?0.其几何意义是：主成分分析相当于对原坐标轴做一次旋转变换，使得新坐标系的第1轴对应于数据变

易的最大方向，第2轴与第1轴正交，且对应于数据变易的第二大方向，依次类推. 3)基本步骤

（1）确定分析变量，收集原始数据；设原始数据矩阵为X?(xij)n?p其中xij表示第i个样品（对象）在第j个变量上的取值。

（2）在进行主成分分析之前，要检验该样本矩阵是否适合于主成分分析.KMO检验是检验变量之间偏相关关系的统计量，用于检验变量间的偏相关系数是否过小. KMO统计量越接近于1，说明各变量间的偏相关系数越大，KMO统计量大于0.9，效果最好；如果统计量小于0.6，则不适合于做主成分分析.Bartlett球形检验是检验相关矩阵是否是单位矩阵，即各变量是否各自独立.

9

（3）对原始数据进行标准化，即令

*xij?xij?xjsj （6）

其中xj,sj分别为第j列元素的样本均值和样本标准差，即

1n1nxj??xij,sj?(xij?xj)2 ?ni?1n?1i?1*则X*?(xij)n?p为标准化的样本资料库.

（4）由标准化后的数据矩阵求协方差矩阵?，或者由原始数据矩阵求相关系数矩阵R.这两种方法结果相等.本文采用直接计算原始数据的相关矩阵的方法（对于数量级差别较大或者有量纲的数据宜适用）.设原始数据X的相关系数矩阵为

?r11r12...r1p???rr...r21222p? （7） R???...????r?r...rnp??n1n2rij(i,j?1,2,?,p)为原变量xi与xj的相关系数，rij?rji，其计算公式为

rij??(xk?1nki?xi)(xkj?xj)n?(xk?1n （8）

ki?xi)2?(xkj?xj)2k?1（5）计算R的特征根和特征向量；

根据特征方程?E?R?0得R的特征根为?i(i?1,2,...,p)，将特征根按照从大到小的顺序排列，排列后的特征根不妨仍然表示为?1??2?...?p?0.同时可得对应的特征向量u1,u2,...,up，将他们标准正交化，u1,u2,...,up称为主轴.

（6）计算所有变量的方差贡献率及累计方差贡献率；?i的方差贡献率为

ei??i??i?1p?i?1,2,?,p? （9）

i?i的累计方差贡献率为

10

Ei???k?1pi?1ik ?i?1,2,?,p,k?1,2,...,m,m?p? （10）

??i（7）确定主成分的数目m. 方法有：①一般取累计贡献率达85%—95%的主成分； ②选用所有?i?1的主成分；③累计特征值乘积大于1的主成分；④画出特征值变化曲线，以转折点位置为标准判断.本文采用累计贡献率达85%—95%的主成分.

（8）确定主成分函数表达式模型. 设m个主成分对应的特征向量分别为

A1、A2、、...Am，其中Aj??a1j主成分yj的函数表达式为

a2j...apj?，akj表示Aj的第k行的元素，则第j个

T?x1???x2?T?yj??Aj??a?...??1j??x???p?a2j?x1???px2??...apj??akjxk （11） ?...??k?1???x??p?（9）提炼主成分yj的抽象意义.由xk与yj的相关系数bkj的大小可以确定yj主要与哪几个变量显著相关，然后根据这几个变量的实际意义提炼yj的抽象意义.

（10）检验主成分模型.根据n个样本的m个主成分的函数值，通过计算m个主成分y1,y2,...,ym的相关系数就可以检验m个主成分是否线性无关.如果两个主成分的相关系数为0，则说明这两个主成分线性无关，模型有效；否则线性相关，模型无效.

（11）求主成分函数值。将各样本标准化数据xk代入（7），可以求得各样本的第

j个主成分yj的函数值.

4)模型求解

（1）收集原始数据矩阵X.

本文选取了某地区的月平均气压的平均值、月最高气压的平均值、月最低气压的平均值、月平均气温的平均值、月最高气温的平均值、月平均气压的平均值8项指标，并分别记为x1，x2??x8. 每个指标有48个数据（见附件1）。

使用SPSS软件进行求解（见附录7）。

（2）将原始数据标准化，（SPSS内部计算）.

（3）求原始数据的相关系数矩阵R，如图7所示.

11

图7 相关系数矩阵

图8因子分析检验图

从图8看出，表格的第一行为检验变量间偏相关程度的KMO统计量，其值在0.6之上才适合做主成分分析，效果显著，如果小于0.6，效果不显著，不适合做主成分分析。下面的三行为球形检验的结果，球形检验原假设的变量是不相关的，显然只有拒绝原假设的情况下数据才适合做因子分析。本例中KMO值为0.720，球形检验显著，两个条件都满足，变量间相关程度大，适合做因子分析。

（4）计算矩阵R的特征根、各因子的方差贡献率及累计方差贡献率，并确定主成分的个数.如图9所示。

图9 R特征值及其累计方差贡献率

从图9中可以看出，第一、第二主成分对方差的累计贡献率达到95.461%，它们分别对应着原样本数据点数据变异的最大、次大方向，是原变量系统的一个最佳整合，从而我们可以以95.461%的精度将变量的有效维数从8维降至2维.因此可以将前2个因子作为主因子.

（5）确定主成分函数表达式模型，因子得分系数矩阵如图10所示.

12

图10 因子得分系数矩阵

设2个主成分分别为y1,y2，则建立模型为

?y1??0.184x1?0.183x2?0.185x3?0.175x4?0.186x5?0.164x6?0.096x7?0.121x8 （12） ??y2?0.048x1?0.042x2?0.057x3?0.020x4?0.055x5?0.015x6?0.536x7?0.567x8 其中x1,x2,??x8.均为原变量经过均值为0，方差为1标准化后的变量. （6）对主成分y1,y2的意义进行解释。图11给出了原变量与第1、第2主成分的相关系数.

图11 旋转后的因子载荷矩阵

第一主成分y1，与原变量x1(平均气压的平均值)、x2（最高气压的平均值）x3（最低气压的平均值）、x6（最低气温的平均值）的相关系数的绝对值都超过了0.948，因此它是一个反映气温和气压的综合因子，我们称之为气压温度因子.

第二主成分y2，与原变量x7（月平均相对湿度的平均值）的相关系数为0.925、

x8（月最低相对湿度的平均值）的相关系数为0.948，其余的都不超过0.266，因此它是一个反映相对湿度的因子，称为湿度因子. （7）计算2个主成分的函数值.

将48个经过标准化的数据xi*代入模型yj，可以得到48个地区的主成分yj的函数

13

值，结果如表1所示.

表1 48个月对应的主成分函数值 序号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 y1 -1.6817 y2 0.5485 z 487 528 784 858 1082 1137 1043 1080 1206 1265 1202 1297 1209 1716 1787 1598 1755 1533 2078 1322 1384 1331 1475 1472 序号 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 y1 -1.121 y2 1.77 z 946 799 739 861 855 788 752 962 828 769 667 623 841 1280 1671 1553 1715 1618 1621 1716 1647 1709 1615 1719 -1.601 -0.061 -0.675 0.0076 0.1144 -0.917 0.9219 -2.098 1.3417 0.0584 1.3521 0.4715 0.9402 1.6644 0.3575 1.5115 0.2595 -0.727 -1.091 1.1667 -1.331 0.3153 -1.453 0.2577 -1.083 0.8107 -0.823 0.1507 -0.401 -0.076 0.7358 -0.659 0.7824 0.7099 1.2455 0.6445 1.3537 -0.133 0.5816 0.9467 -0.246 0.4312 -0.531 -0.633 -0.593 -1.472 -0.7683 -0.2082 -0.4507 -0.2166 0.152 0.9259 1.5239 0.5675 -0.1318 -1.1893 -1.5909 0.7554 0.2097 1.0055 0.192 0.5304 1.1598 -1.7758 1.5284 -0.4272 -0.6524 -0.9464 -1.5177 -0.4902 -1.3599 -1.3556 -0.1654 -2.0187 0.1845 0.6764 1.1459 0.5359 -0.1899 -0.9077 -1.03 1.8659 0.7356 1.1416 0.9149 0.5766 0.9955 -1.3149 1.6182 -0.2525 -1.036 -0.9904

（8）检验主成分模型.

由于主成分分析的4个条件中的前3个（每一个主成分都是各原始变量的线性组合；主成分的数目大大少于原始变量的数目；主成分保留了原始变量的绝大多数信息），只要检验4个主成分是否相关即可.

由步骤（7）计算的2个主成分的得分矩阵Y??yij?48?2，求矩阵Y的协方差矩阵如图12所示.

14

图12 因子得分的协方差矩阵

从图12可以看出，主成分得分的协方差矩阵为单位矩阵，说明提取的2个主成分是互不相关的.满足假设的条件，模型和结果有效. 5.2.2多元非线性回归分析[4]

将48个月的发病率作为因变量，记作z，将发病率的48个数据填入表1中。下面寻找发病率z与主成分y1,y2的关系式，这需要使用多元非线性回归分析方法。

经过反复试验探索，找到的非线性回归模型为

??aya4y2az?exp?a1y1?a2y2?31??5?a6? （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的系数为负值可得，气压温度因子与月平均气压、月平均最高气压、月平均最低气压具有负相关性；

（2）由x4,x5,x6的系数为正值可得，气压温度因子与月平均温度、月平均最高温度、月平均最低温度具有正相关性；

（3）由x7,x8的系数为负值可得，气压温度因子与月平均相对湿度、月平均最低相对湿度具有负相关性； 由第2个方程可知：

（4）由x1,x2,x3的系数为正值可得，湿度因子与月平均气压、月平均最高气压、月平均最低气压具有正相关性；

（5）由x4?0,x5?0,x6?0可得，湿度因子与月平均温度、月平均最高温度成负相关性，与月平均最低温度具有正相关性；

（6）由x7,x8的系数为正值可得，湿度因子与月平均相对湿度、月平均最低相对湿度具有正相关性；

15

5.3高危人群预警分析（问题3）

首先预测2011年的气象状况，然后预测未来2011年的高危人群的发病率和发病时间，最后提出预警和干预措施。 5.3.1未来2011年气象预测

1）气压、温度状况预测

观察月平均气压x1、月平均最高气压x2、月平均最低气压x3、月平均温度x4、月平均最高温度x5，月平均最低温度x6的历史数据，发现随时间做周期性变化，于是建立余弦函数模型x?Acos(?t??)?B，利用过去36个月的历史数据进行参数估计，然后使用2009年的12个数据进行预测，评估误差并检验模型的可靠性。最后预测出2010年的12个月的数据。建模结果见表2.（MATLAB程序见附录8）。2011年12个月的预测结果见表3.

表2 模型参数估计结果 变量 表达式 参数拟合值 平均相对误差 0.15% x1 x?Acos(?t??)?B A??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.8618x2 x?Acos(?t??)?B 0.14% x3 x?Acos(?t??)?B 0.15% x4 x?Acos(?t??)?B 11.77% x5 A??12.0959,??0.5198, x?Acos(?t??)?B ???6.8425,B?21.0561x?Acos(?t??)?B 7.75% x6

A??12.2819,??0.5224, ???0.6683,B?13.548234.71%

16

表3 2011年12个月的气压、温度预测值 月份 平均温度 平均压力 最高压力 最低压力 最高温度 最低温度 月份 平均温度 平均压力 最高压力 最低压力 最高温度 最低温度 1 4.9978 1026.1 1018.2 1014.7 23.7426 16.0265 7 28.6862 1004.9 1017.6 1011.7 18.0985 10.9839 2 5.2969 1024.3 1023..7 1019.5 17.5302 9.6937 8 28.4803 1006.7 1012.1 1006.9 24.3144 17.3191 3 8.6662 1020.1 1027.7 1022.6 12.249 4.3892 9 25.1901 1010.9 1008.1 1003.8 29.6697 22.6484 4 14.2111 1014.7 1029 1023.2 9.294 1.5279 10 19.689 1016.2 1006.6 1003.2 32.75 25.5502 5 20.4596 1009.6 1027.4 1021.1 9.4456 1.8731 11 13.4374 1021.4 1008.2 1005.2 32.7417 25.2503 6 25.7531 1006 1023.2 1016.9 12.6638 5.3328 12 8.0948 1025.1 1012.3 1009.4 29.647 21.8289

2）湿度状况预测

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

（1）自相关系数

原始序列X(0)=?X(0)(1),X(0)(2),...,X(0)(n?)的各阶自相关系数反映已知数据对未

来数据的影响程度. 各阶自相关系数为

n?wrw????Xk?1(0)(0)(0)??(k)?X(0)?X(k?w)?X??????Xk?1n(0)(k)?X(0)2 （13）

??式中，

X对各阶自相关系数归一化得，

(0)1n(0)=?X(k) （14） nk?1?w=rw?rw?1t,w?1,2,...,t （15）

w17

?w可作为各阶步长的马尔柯夫链权重，t是按预测需要计算的最大阶数，一般取

rw?0.3.根据rw?0.3可以确定转移步数w.

（2）加权马尔柯夫模型

①状态划分。设划分的m个湿度区间为

E1??a11,a21?,Ei??a1i,a2i?,i?2,3,...,m?1 Em??a1m,a2m?其中，a11尽可能小，a2m尽可能大.，如果

?(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)(w)w)?p11?p12...p1(m?(w)(w)(w)?pp...p222m?R(w)??21 （16） ?...??(w)(w)(w)?pp...pm2mm??m1已知每一步的概率转移矩阵和每一步的初始状态，则马尔柯夫链就可以确定. ③预测值计算

选取距离预测年最近的t(t?m)个年份，按照距离预测年由近到远，转移步数w分别为1,2,...,t，以这几年的相对误差所对应的状态为初始状态，不妨设第1,2,...,t年所

(1)(2)对应的初始状态分别为E1',E2',...,Et('t)，其中，w'??1,2,...,m?.例如，当2'?5时，(2)E2'=E5(2)，说明距离预测年第2年的状态是第5状态.在转移步数w对应的转移矩阵(w)(w)(w)(w)p(w)=?pw...pwR(w)中，取起始状态Ew'1，pw'2，，'m?，从而组成新的概'所对应的行向量w'18

率矩阵

(1)?p1'1?(2)pR??2'1?...?(t)?pt'1(1)p1'2...(2)p2'2...(1)p1'm?(2)?p2'm? （17） ?)?pt('tm?t)pt('2...将矩阵R加权得

(1)(1)(1)??1p1'1?1p1'2...?1p1'm??(2)(2)(2)??p?p...?p22'222'm?R???22'1 （18） ?...??(t)(t)(t)??p?p...?ptt'2tt'm??tt'1将矩阵R?按列求和得

tt?t(w)(w)(w)?p??p1,p2,...,pm?=???wpw'1,??wpw'2,...,??wpw'm? （19）

w?1w?1?w?1?找出向量p的最大分量得

pM?max?p1,p2,...,pm?,M??1,2,...,m? （20）

分量pM所对应的状态EM就是预测年的状态，则该年度的预测值为

?M?a1M?a2M （21） 2（3）计算过程和结果

以预测2010年第1月的数据为例。利用2007、2008、2009年的36个月的历史数据进行建模。

①自相关系数。以根据（13）、（14）、（15）式计算可得各阶的自相关系数，确定最大滞后阶数w?2.各阶自相关系数及权重见表4.

表4 自相关系数及权重

w 1 2

rw

0.4737 0.6625

0.2414 0.3375

?w

②划分的6种状态区间，见表4.

表4 各个状态区间 状态编号 状态区间 E1 [0,60) E2 [60,65) E3 [65,70） E4 [70,75) E5 [75,80） E6 [80,100] ③构造转移概率矩阵

如果有的状态不能从统计表中得到转移概率，则假定它未来转移到各个状态的概

19

率都相等，即都等于表6.

1.根据（15）可得1步和2步内的转移概率矩阵分别见表5和m表5 1步转移概率矩阵R(1)

2/3 0 1/8 0 0 0 1/3 1/5 1/4 1/7 0 0 0 2/5 1/8 3/7 1/9 1/3 0 1/5 0 1/7 4/9 1/3 0 1/5 3/8 2/7 1/3 0 0 0 1/8 0 1/9 1/3 表6 2步转移概率矩阵R(2)

1/3 1/5 1/8 0 0 0 1/3 0 1/8 2/7 1/8 0 1/3 0 0 4/7 1/2 1/3 0 1/5 5/8 0 0 0 0 2/5 1/8 2/7 1/4 2/3 0 1/5 0 1/7 1/8 0

④组成预测年份的新转移概率矩阵.选择离预测年最近的2个年份，转移步数分别为w?1,2，根据（17）式得预测年的转移概率矩阵R37，见表7.

表7 月份 状态 步长 36 35 2 3 1 2 权重 1 0.6625 0.3375 0 0 0 预测年的转移概率矩阵R37 2 1/7 1/8 3 3/7 1/2 4 2/7 0 5 1/7 1/4 6 0 1/8 概率来源 R(1) R(2) 加权求和 0.1368 0.4527 0.0946 0.2737 0.0422 ⑤确定预测年份的状态.预测年的状态向量的最大分量值为0.4527，对应的状态为第3状态，即第37个月的湿度将处于第3状态，湿度67.5.

将第37个月的湿度值67.5放入序列中，同时去掉第1个月的湿度数据，重新构建马尔柯夫链，得第38个月的湿度。以此类推，可得2010年12个月的湿度值，见表8.平均绝对相对误差为7.97%，可靠性高.

表8 2010年湿度预测 月份 实际x7 1 67.5 2 77.5 3 72.5 4 67.5 5 77.5 6 72.5 7 67.5 8 77.5 9 72.5 10 67.5 11 12 77.5 72．5 20

建立M文件夹： %发病人群信息

A=xlsread('c题数据1.xls','sheet1','A2:F58926');X925行-6列 a=size(A);

B=zeros(4,12); for i=1:a(1)

if A(i,4)==2007 if A(i,5)==1

B(1,1)=B(1,1)+1; elseif A(i,5)==2 B(1,2)=B(1,2)+1; elseif A(i,5)==3 B(1,3)=B(1,3)+1; elseif A(i,5)==4 B(1,4)=B(1,4)+1; elseif A(i,5)==5 B(1,5)=B(1,5)+1; elseif A(i,5)==6 B(1,6)=B(1,6)+1; elseif A(i,5)==7 B(1,7)=B(1,7)+1; elseif A(i,5)==8 B(1,8)=B(1,8)+1; elseif A(i,5)==9 B(1,9)=B(1,9)+1; elseif A(i,5)==10 B(1,10)=B(1,10)+1; elseif A(i,5)==11 B(1,11)=B(1,11)+1; elseif A(i,5)==12 B(1,12)=B(1,12)+1;

end 07年各自月份病发者指标 elseif A(i,4)==2008 if A(i,5)==1 B(2,1)=B(2,1)+1; elseif A(i,5)==2 B(2,2)=B(2,2)+1; elseif A(i,5)==3 B(2,3)=B(2,3)+1; elseif A(i,5)==4 B(2,4)=B(2,4)+1; elseif A(i,5)==5 B(2,5)=B(2,5)+1; elseif A(i,5)==6

26

B(2,6)=B(2,6)+1; elseif A(i,5)==7 B(2,7)=B(2,7)+1; elseif A(i,5)==8 B(2,8)=B(2,8)+1; elseif A(i,5)==9 B(2,9)=B(2,9)+1; elseif A(i,5)==10 B(2,10)=B(2,10)+1; elseif A(i,5)==11 B(2,11)=B(2,11)+1; elseif A(i,5)==12 B(2,12)=B(2,12)+1;

end 08年各自月份病发者指标 elseif A(i,4)==2009 if A(i,5)==1 B(3,1)=B(3,1)+1; elseif A(i,5)==2 B(3,2)=B(3,2)+1; elseif A(i,5)==3 B(3,3)=B(3,3)+1; elseif A(i,5)==4 B(3,4)=B(3,4)+1; elseif A(i,5)==5 B(3,5)=B(3,5)+1; elseif A(i,5)==6 B(3,6)=B(3,6)+1; elseif A(i,5)==7 B(3,7)=B(3,7)+1; elseif A(i,5)==8 B(3,8)=B(3,8)+1; elseif A(i,5)==9 B(3,9)=B(3,9)+1; elseif A(i,5)==10 B(3,10)=B(3,10)+1; elseif A(i,5)==11 B(3,11)=B(3,11)+1; elseif A(i,5)==12 B(3,12)=B(3,12)+1;

end 09年各自月份病发者指标 else A(i,4)==2010 if A(i,5)==1 B(4,1)=B(4,1)+1; elseif A(i,5)==2

27

B(4,2)=B(4,2)+1; elseif A(i,5)==3 B(4,3)=B(4,3)+1; elseif A(i,5)==4 B(4,4)=B(4,4)+1; elseif A(i,5)==5 B(4,5)=B(4,5)+1; elseif A(i,5)==6 B(4,6)=B(4,6)+1; elseif A(i,5)==7 B(4,7)=B(4,7)+1; elseif A(i,5)==8 B(4,8)=B(4,8)+1; elseif A(i,5)==9 B(4,9)=B(4,9)+1; elseif A(i,5)==10 B(4,10)=B(4,10)+1; elseif A(i,5)==11 B(4,11)=B(4,11)+1; elseif A(i,5)==12 B(4,12)=B(4,12)+1;

end 10年各自月份病发者指标 end end

C=mean(B)/mean(mean(B)),%季节指数 附录7

主成分分析：

设置步骤如下： （1）打开因子分析.

28

（2）将要分析的变量移入.

（3）选中“初始解”、“相关系数”和“KMO检验和Bartlett球形检验”.

（4）选中“主成分分析方法”、“相关系数矩阵”、“未经过旋转的因子解（因子载荷矩阵）”、“碎石图”、“变量提取数目”.

（5）（如果需要旋转的话）选中“方差最大化旋转方式”、“旋转解”、“因子载荷图”.

点击“OK”就可以输出计算结果.

附录8 程序： x=1:48;

29

y=[1.87 5.01 7.95 11.29 18.49 22.12 26.29 26.61 21.66 16.05 8.93 5.2 0.83 -0.59 7.01 12.1 17.04 21.05 26.81 24.94 22.12 16.82 8.65 2.71 -0.37 5.71 6.17 11.66 16.72 22.5 25.1 25.26 21.72 16.5 7.97 2.54 1.03 3.81 5.44 9.02 17.13 20.76 25.84 26.97 22.67 14.67 8.69 3.08 ];

f=inline('a(1)*cos(a(2)*x+a(3))+a(4)','a','x'); a=lsqcurvefit(f,[100;pi/6;-pi/4;30],x,y) x1=1:0.01:48;

y1=a(1)*cos(a(2)*x1+a(3))+a(4); plot(x,y,'o',x1,y1,'r')

y2=a(1)*cos(a(2)*x+a(3))+a(4); y3=(y2-y)./y;%输出相对误差 y4=abs(y3);%输出绝对相对误差

y5=mean(y4),%输出平均绝对相对误差 y6=max(y4);%输出最大绝对相对误差 t=49:1:60;

y7=a(1)*sin(a(2)*t+a(3))+a(4)

30

本文来源：https://www.bwwdw.com/article/f0da.html