数据分析与处理

一个地区旱涝灾害受降水量的影响，同时对旱涝灾害的估计也是预备救援物资的依据。表1是某地区45年来的降水量的数据（单位：毫米），请依次完成以下问题：

（1） 根据表1中的降雨量数据，给出第46年的降水量预测值，并说明结果的置信度。

（2） 根据表1中的降雨量数据，建立数学模型将旱涝灾害分级，并说明第46年的灾害分级结果； （3） 建立数学模型说明在较长时间内，该地区准备抗洪和抗旱物资的策略。

表1 某地区45年来降水量数据（单位：毫米）

年份 1月 2月 3月 4月 5月 6月 7月 8月 9月 10月 11月 12月 年总计 1 2 3 5 6 7 8 9 10 11 12 13 14 15 16 17 10 20 24 14 43 52 63 23 19 16 0 35 18 18 52 82 55 92 230 127 341 104 127 93 162 106 46 160 201 48 130 140 90 165 43 78 55 160 72 85 2 28 5 15 58 44 63 77 29 0 22 30 155 68 24 56 32 10 145 61 86 68 46 21 46 33 1 4 142 13 102 66 100 71 75 6 57 50 160 26 1 35 20 4 31 32 60 0 26 39 26 11 38 51 2 1284 1030 1090 2031 1193 975 1322 1410 19 1017 1034 1621 1095 1344 8 704 1158 77 169 47 114 94 4 103 25 224 351 470 568 53 159 103 111 4 14 185 107 245 196 42 48 98 269 181 60 148 153 242 295 112 87 243 175 519 79 247 46 32 0 17 115 306 113 23 104 30 1 32 71 49 19 31 160 66 216 24 88 86 28 67 110 199 119 26 160 105 13 136 148 203 385 416 108 77 223 126 40 183 9 166 327 60 81 87 32 73 0 91 196 230 379 44 168 44 100 111 110 65 139 127 155 76 101 68 145 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 67 10 41 28 28 63 10 26 35 28 42 51 28 37 34 8 10 58 15 0 120 67 53 84 138 82 200 57 80 27 29 73 83 19 74 6 47 43 61 23 1 26 101 204 23 56 124 101 45 63 0 83 169 14 408 117 101 227 26 138 36 4 9 32 38 24 32 132 1 44 22 20 44 48 6 20 91 130 35 66 43 45 84 0 111 78 2 14 104 2 35 8 26 1 37 59 10 46 13 42 6 0 3 18 81 18 41 0 5 22 39 45 33 963 1727 1208 776 10 1215 943 1308 871 1178 7 976 1601 1132 1596 1873 11 1005 1030 1424 1302 1629 1306 1782 1094 44 106 162 78 65 205 530 480 43 72 48 41 20 18 14 33 86 144 206 251 153 84 126 133 75 13 193 59 101 107 136 50 51 90 102 143 2 141 173 50 136 194 101 138 72 86 98 100 269 77 204 162 26 27 92 32 217 65 220 137 250 101 186 19 119 19 18 129 253 193 117 112 121 81 123 119 386 34 219 38 185 73 127 15 19 47 10 92 56 146 275 286 421 80 168 34 113 112 68 143 413 201 258 138 74 451 122 77 138 62 165 222 91 32 45 95 1 30 168 191 385 321 114 151 52 100 97 74 149 50 80 253 23 128 52 104 148 195 163 173 219 36 302 225 81 310 160 38 36 62 67 104 47 180 18 88 148 180 355 134 163 119 83 197 145 217 129 117 43 114 125 171 211 193 716 111 27 223 102 143 331 43 44 45 98 23 79 88 127 123 255 135 190 119 201 90 43 63 103 80 92 316 37 118 0 41 203 257 223 169 153 58 25 104 129 46 0 28 27 2 1551 1020 1274

【问题1】：

 插值与拟合（要求降雨量与年份有确定的联系，满足吗？）  回归分析方法（降雨量与年份有因果联系吗？）  时间序列方法（可行吗？）

 灰色预测（对数据进行累加后，可能出现降雨量总合与时间的联系？） 使用灰色步骤： 【Step 1】：数据检验与预处理：级比检验和轻易变换 【Step 2】：建立GM(1,1)模型并求解 【Step 3】：模型检验，确定可行性 【Step 4】：预测预报 【问题2】：

【Step 1】：建立旱涝灾害分类模型（聚类模型）； 【Step 2】：求解模型，从而将它们分类； 【Step 3】：建立判别分析模型，将预测结果归类。

两种选择：第一选择：只对年降雨量分类：特干(1)、较干、正常、较涝、特涝（可以采用欧氏做聚类分析，分成五类）

Y=pdist(A(:,13)); Z=linkage(Y,'ward'); [H,T]=dendrogram(Z); c=cophenet(Z,Y);

T=cluster(Z,'maxclust',5);

第二选择：对所有数据做聚类分析，分成合适的类（较复杂）。 分类结果： 年份 一月 10 20 24 103 14 43 52 63 31 23 19 16 0 35 18 18 52 43 67 10 41 二月 82 77 47 94 53 14 42 48 160 66 26 24 88 86 28 65 0 44 67 53 三月 55 169 114 25 159 185 60 98 87 216 160 13 77 91 40 44 139 120 106 86 78 四月 92 93 46 224 103 107 148 269 243 67 105 136 223 196 183 100 127 84 162 144 84 五月 230 162 160 351 111 245 153 181 175 110 79 148 126 230 44 111 155 138 65 206 126 六月 127 106 201 470 4 196 242 17 519 199 247 203 9 379 168 110 76 82 205 251 133 七月 341 48 90 568 55 72 295 115 46 119 32 385 166 60 81 87 101 200 530 153 43 八月 104 130 165 43 160 78 112 306 0 23 30 416 327 32 73 0 68 57 480 13 48 九月 127 140 85 2 28 5 15 113 58 104 1 108 32 71 49 19 145 80 27 193 29 十月 44 63 77 29 0 22 30 155 68 24 56 32 10 145 61 86 68 23 1 26 101 十一月 46 21 46 33 1 4 142 13 102 66 100 71 75 6 57 50 160 32 38 24 32 十二月 26 1 35 20 4 31 32 60 0 26 39 26 11 38 51 2 104 2 35 8 年总量 1284 1030 1090 2031 1193 975 1322 1410 19 1017 1034 1621 1095 1344 8 704 1158 963 1727 1208 776 分类结果 2 5 1 4 1 5 2 2 3 5 5 3 1 2 5 5 1 5 4 1 5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 28 28 63 10 19 26 35 28 42 51 28 37 34 8 10 58 15 67 47 18 98 23 79 59 90 50 51 119 18 19 38 34 73 15 19 47 10 52 92 104 180 114 27 88 90 43 101 102 50 65 72 129 86 81 185 113 127 30 52 97 74 104 50 88 83 125 223 127 63 41 107 143 136 220 98 253 77 123 56 112 68 168 100 80 149 148 36 148 197 171 102 123 103 203 136 2 194 137 100 193 204 119 146 34 143 191 74 253 62 195 302 180 145 211 143 255 80 257 75 141 101 250 269 117 162 386 275 219 413 385 451 77 165 163 225 355 217 193 331 135 92 223 72 173 138 101 26 112 27 92 286 80 201 321 122 138 222 173 81 134 129 716 190 316 169 41 32 20 186 18 121 14 33 421 168 258 114 91 32 23 219 310 163 117 111 43 119 37 153 73 217 83 19 74 6 47 43 61 138 151 45 95 128 1 160 119 38 36 62 201 118 0 204 23 56 124 101 45 63 0 83 169 14 408 117 101 227 26 138 36 4 9 58 25 104 132 1 44 22 20 44 48 6 20 91 130 35 66 43 45 84 0 111 78 2 14 129 46 0 26 1 37 59 10 46 13 42 6 0 3 18 81 18 41 0 5 22 39 45 33 28 27 2 10 1215 943 1308 871 1178 7 976 1601 1132 1596 1873 11 1005 1030 1424 1302 1629 1306 1782 1094 1551 1020 1274 5 1 5 2 5 1 5 5 3 1 3 4 1 5 5 2 2 3 2 4 1 3 5 2 【问题 3】：

【Step 1】：建立平稳马尔可夫链模型

根据状态转移频率建立转移矩阵，建立马尔可夫过程：

WWPw1w2w3w4w51 w0i【Step 2】：求解该模型：计算特征向量。 【Step 3】：解读数据并确定整体的物资准备策略： 【Step 4】：利用问题1的模型连续计算几年的结果，并判别特殊状态。

第一步：探索性数据分析（统计方法）

1．做出相关图形：散点图、直方图、频数表

% 1. 绘制：散点图、直方图、盒子图、经验累加分布函数图

plot(A(:,1),'*') %绘制散点图 hist(A(:,1)) %绘制直方图 boxplot(A(:,1:12),1,'+',0)%盒子图 cdfplot(A(:,1)) %经验累加分布函数图

经过观察无明显规律。

2．计算相关的统计参数：频数表、均值、方差、标准差、偏度、峰度、协方差矩阵等

tabulate(A(:,1)) % 频数表 M=mean(A); % 均值 V=var(A); % 方差

S=std(A); % 标准差 y=skewness(A) % 偏度

【解读】：偏度用于衡量样本均值的对称性，若偏度为负，则数据均值左侧的离散型比右侧的强；若偏度为正，则右侧的离散性比左侧的强。严格对称分布的偏度为0。 【结果】：0.99934 全部右偏。 k=kurtosis(A) % 峰度

【解读】：峰度是分布形状的另一种度量，若比3（正态分布的峰度）大得多，表示分布有沉重的尾巴，说明样本中有较多远离均值的数据。 【结果】：3.6784

1.0284 0.70601 0.56734 0.403 0.6779 1.9055 1.3992 0.558 2.3232 0.82274 1.1327 0.62368。

4.3238 3.1295 2.5877 3.0065 2.7328 6.5922 4.2458 2.4431 10.269 2.8131 4.3174 2.8452；

10月分布最分散。

[R,P]=corrcoef(A(:,1:12)) % 相关系数

【解读】：R返回相关系数矩阵；P用于检验假设：没有相关性(0.05)，如果P(i,j)较小，小于0.05，则相关性R(i,j)显著。 【结果】：

相关系数表格

1 0.149 0.27593 0.35991 0.087625 0.36359 0.10318 -0.19965 0.1119 0.17558 0.149 1 0.18386 0.23931 0.1944 0.085131 -0.10886 0.0692 -0.13387 0.1292 -0.070286 -0.1383 1 -0.31307 -0.12404 -0.024753 -0.21859 -0.00774 0.093672 -0.23614 0.27593 0.18386 -0.31307 1 0.13812 0.05388 0.18749 0.06594 -0.156 0.14996 0.35991 0.23931 -0.12404 0.13812 1 0.034333 0.20326 0.14232 0.014662 0.087625 0.1944 0.05388 0.034333 1 0.045256 -0.1498 0.36359 0.085131 0.18749 0.20326 0.045256 1 0.10318 -0.10886 0.06594 -0.060145 -0.058771 0.32935 1 -0.0293 -0.19965 0.0692 -0.156 0.14232 -0.1498 0.1119 -0.13387 -0.23614 0.14996 0.014662 0.17558 0.1292 -0.026847 -0.23735 -0.077931 -0.011184 0.15171 0.1106 1 0.050781 -0.062304 0.29384 -0.018665 0.18425 0.2104 -0.24315 -0.14561 -0.18325 -0.026165 -0.070286 -0.1383 -0.024753 -0.21859 -0.00774 0.093672 -0.00043908 -0.32115 0.0018883 -0.061987 -0.0634 -0.1313 -0.037819 -0.0293 1 0.059968 0.15171 0.059968 1 0.1106 -0.060145 -0.058771 0.32935 0.0018883 -0.1313 -0.0634 -0.037819 -0.00043908 -0.026847 -0.23735 -0.061987 -0.077931 -0.011184 0.050781

-0.062304 -0.32115 0.29384 -0.018665 0.18425 0.2104 -0.24315 -0.14561 -0.18325 -0.026165 1 相关性检验表

1 0.32687 0.638 0.0669 0.015166 0.56707 0.014084 0.50001 0.18855 0.428 0.24862 0.74043

【解释】：

cov(A(:,1:12)) % 协方差阵

[tab,chi2,p]=crosstab(A(:,11),A(:,12));% 列联表检验性

0.32687 1 0.393 0.22666 0.11336 0.20068 0.5782 0.47656 0.55793 0.38063 0.397 0.68432 0.638 0.393 1 0.036262 0.4169 0.87178 0.14913 0.96632 0.051 0.11837 0.99772 0.031473 0.0669 0.22666 0.036262 1 0.36557 0.7252 0.21746 0.66693 0.30338 0.329 0.86103 0.050091 0.015166 0.11336 0.4169 0.36557 1 0.82284 0.18052 0.69472 0.35102 0.92384 0.113 0.90314 0.56707 0.20068 0.87178 0.7252 0.82284 1 0.76784 0.70136 0.32603 0.99018 0.68584 0.22567 0.014084 0.5782 0.14913 0.21746 0.18052 0.76784 1 0.027161 0.6788 0.393 0.61087 0.16537 0.50001 0.47656 0.96632 0.66693 0.69472 0.70136 0.027161 1 0.80518 0.84848 0.94187 0.10752 0.18855 0.55793 0.051 0.30338 0.35102 0.32603 0.6788 0.80518 1 0.69557 0.31982 0.339 0.428 0.38063 0.11837 0.329 0.92384 0.99018 0.393 0.84848 0.69557 1 0.46952 0.22822 0.24862 0.397 0.99772 0.86103 0.113 0.68584 0.61087 0.94187 0.31982 0.46952 1 0.853 0.74043 0.68432 0.031473 0.050091 0.90314 0.22567 0.16537 0.10752 0.339 0.22822 0.853 1 【解读】：chi2为统计量，用于检验表中行和列的性。标量p为检验的显著性水平。当p接近于0时，可以拒绝零假设，认为行和列之间是不的。

【结果】：p=0.4634。11月和12月降雨量是各自的。

3．统计推断：判定数据服从怎样的分布，并作分布检验

【变量分布形态的估计】：

 频数分布表和频率直方图

 频数表：tabulate(A(:,1))

 直方图：hist(A(:,1))

 带正态密度曲线的直方图：histfit(A(:,1))  经验分布函数

 cdfplot(A(:,1))

 五数概括与Box图

 Boxplot(A(:,1))

【变量分布参数的估计】：

 矩估计：原点矩、方差  极大似然估计  区间估计

【变量分布形态的检验】：

 假设检验的步骤

 提出原假设H0和备择假设H1；

 选取一个适当的统计量T，并写出相应的检验准则；  给定显著性水平，并求出H0的拒绝域W；

 由样本算出检验统计量T的实测值，判断其是否落入拒绝域。  K.Pearson-Fisher检验



2拟合优度检验：关于变量X分布形态的某种先验知识或猜测是否为真的统计推断方法。

 用极大似然法估计分布参数：1,2,,r；  对假设H0:FX(x)F0(x;1,2,,r)进行检验。

【例】：对一月份的降雨量数据（数据表中第一列）做正态性检验。

分析：这是一个正态拟合检验问题。检验的原假设为H0:X1N(,2)，其中参数,2未知。 第一步：进行未知参数的极大似然估计

[mu,sigma]=normfit(A(:,1));%输出：mu=35.7333;sigma=23.4370。

于是原假设修正为：H0:X1N(35.7333,23.4372)。 第二步：样本数据分组

[f,med]=hist(A(:,1)); f_med=[f',med'] 运行结果：

6 9 8 6 5 4 4 1 0 2 5.15 15.45 25.75 36.05 46.35 56.65 66.95 77.25 87.55 97.85 利用hist指令自动分为10组，并统计各组频数。由计算结果知道，前后三组数据的频数偏小，故将后三组数据进行合并，这样可得8组数据。这8组数据所属的数据组的区间边界值如下：

a=[]; for k=1:7

aa=(med(k)+med(k+1))/2; a=[a,aa]; end

a=[-inf,a,inf]

输出如下：

-Inf 10.3000 20.6000 30.9000 41.2000 51.5000 61.8000 72.1000 Inf 第三步：统计经验频数

经验频数在第二步中已经给出，只需将最后三组合并即可。 f=[f(1:7),f(8)+f(9)+f(10)];

输出结果为：6 9 8 6 5 4 4 3

第四步：计算理论频数

pest=[]; for i=1:8

pp=normcdf(a(i+1),mu,sigma)-normcdf(a(i),mu,sigma); pest=[pest,pp]; end thef=45*pest

输出结果：6.2514 5.4142 7.1582 7.8259 7.0750 5.21 3.2695 2.7166。 第五步：计算检验统计量的观测值

chi2est=sum((f-thef).^2./thef) 输出结果： 4.0256。 第六步：检验决策

k=8;%分组数量 r=2;%参数数量 alpha=0.20;%检验水平 df=k-r-1;

refcr=chi2inv(1-alpha,df);%计算拒绝域临界值 p=1-chi2cdf(chi2est,df);%检验的p值 if chi2est>refcr h=1; % 拒绝原假设 else

h=0; % 接受原假设 end

alpha,h,p,chi2est,refcr

输出结果：0.10，0，0.57，4.0256，9.23。

计算结果表明，在0.10显著水平下，h=0保留原假设H0，即拟合优度检验认为一月份的降雨量X12N(35.7333,23.4572)。由最小

显著性概率p=0.57表面，当前样本数据下不能拒绝原假设。

 列联表的性检验7

[tab,chi2,p]=crosstab(A(:,11),A(:,12));

输出结果：p=0.4634。11月和12月降雨量是各自的。

 Kolmogorov-Simirnov 检验

原理：参见概率统计教材中的假设检验。 函数：[h,p,stats,cv]=kstest(x,cdf,alpha,tail) 参数说明：x：样本数据向量；

cdf：检验的原假设所指定的分布形式（具体引用为变量的累计分布函数，缺省时cdf=[]，表示拟合标准正态分布）； alpha：检验的显著性水平（缺省时为0.05）； tail：备择假设类型的标示值； h：检验决策；

p：拒绝原假设的最小显著概率； stats：检验统计量的值； cv：拒绝域的临界值。

x=(A(:,1)-mu)/sigma;

[h,p,stats,cv]=kstest(x,[],0.10,0);

输出结果：h=0，p=0.3142，stats=0.1404，cv=0.1718。 结果表面，接受原假设，即服从正态分布。

 正态性检验（相关原理参考有关教材）

 概率纸法：normplot(x)

 Lilliefors法：[h,p,stats,cv]=lillietest(x,alpha,tail)  Jarque-Bera法：[h,p,stats,cv]=jbtest(x,alpha,tail)

4．可以用于回答问题吗？

【问题1】：预测？不可以！需要新的预测方法 【问题2】：分类模型 【问题3】：用处？