缺失数据地补充及异常数据地修正.docx

1、缺失数据地补充及异常数据地修正数模培训作业论文题目 缺失数据地补充及异常数据地修正 缺失数据地补充及异常数据地修正摘 要数据处理贯穿于社会生产和社会生活的各个领域。数据处理技术的发展及其应用的广度和深度，极大地影响着人类社会发展的进程。数据补充，异常数据的鉴别及修正，在各个领域也起到了重要作用。 针对第一问，我们采用了两种模型。第一种是一元多项式回归模型，适用于只有一种自变量的情况。利用我们找到的数据，首先作出散点图，观察其形状，决定拟合多项式的次数，得出拟合曲线与拟合多项式。之后算出均方根误差验证拟合效果，均方根误差较小，说明拟合曲线与源数据吻合得较好。若x1=37.25,x2=41.75,

2、x3=44.5时,y的数据缺失,将x1、x2、x3的值带入拟合多项式，算出缺失值y1=3.3257,y2=2.0437,y3=4.6002，即可补充缺失数据。第二种是多元线性回归模型，适用于有多个自变量的情况。利用我们找到的数据，首先作出散点图，之后作多元回归，求出多元线性回归多项式，以及置信区间。作出残差分析图验证拟合效果，残差较小，说明回归多项式与源数据吻合得较好。若x1=0.055，x2=0.025时，y的数据缺失，则将x1，x2带入回归多项式，算出缺失值y=0.052792。类似地，若x1=0.110,x2=0.045时，y的数据缺失，则将x1，x2带入回归多项式，算出缺失值y=0.0

3、70212，即可补充缺失数据。 针对第二问，我们使用了异常值检验中标准差未知的t检验法。首先绘制火柴棒图观察可疑测定值，可得到可疑值为第6,9,13,23,26,29,35,36,45,53行的数据。将除可疑测定值以外的其余测定值当做一个总体，并假设该总体服从正态分布。由这些测定值计算平均值与标准差s，而将可疑值分别当做一个样本容量为1的特殊总体。如果可疑值与其余测定值同属于一个总体，则它与其余测定值之间不应有显著性差异。检测统计量为：，假设可由标准差s替代来进行检验,则检测统计量可视为：。若统计量值大于相应置信度下的t检验法的临界值（该临界值通过查表法得出），则将可疑值判为异常值。通过计算我

4、们发现，上述可疑值都是异常值。针对第三问，我们采用了分段线性插值、三次样条函数插值以及分段三次Hermite插值法来修正数据异常。同时也需利用外插法修正最后一个数据的异常。之后利用第二问中的t检验判断修正后的值是否仍为异常值。检验结果显示：分段线性插值、三次样条函数插值所修正的第7个和第9个数据仍为异常值，而分段三次Hermite插值所得到的修正值全都不是异常值。所以运用分段三次Hermite插值得到的结果较准确。关键词：一元多项式回归、多元线性回归、t检验法、分段线性插值法、三次样条函数插值、分段三次Hermite插值一、问题重述1、关于数据缺失时应该怎样地把缺失数据补充上来？2、数据完整，

5、但是数据出现异常，如何给出模型找出异常数据？3、异常数据如何修正？2、模型假设1、假设只有因变量存在数据缺失，而自变量不存在缺失。2、利用t检验法时，将除可疑测定值以外的其余测定值当做一个总体，并假设该总体服从正态分布。3、假设可由样本值计算标准差s替代来进行检验。3、符号说明符号含义RMSE均方根误差样本平均值s样本标准差可疑测定值总体标准差k检测统计量置信度置信度下的t检验法临界值n样本容量4、对问题一的分析和处理4.1 一元多项式回归模型1当有缺失的一组数据只有一个自变量时，可以考虑使用一元多项式回归模型。我们采用的数据如下表：表1首先使用matlab作出散点图，观察其形状，决定拟合多项

6、式的次数，运用matlab编程（matlab程序见附录一）得出拟合曲线与拟合多项式为： z=0.16599*x2-13.387*x+271.62。 拟合曲线如下图所示：图1 拟合曲线之后算出均方根误差RMSE验证拟合效果，均方根误差RMSE=0.13931较小，说明拟合曲线与源数据吻合得较好。若x1=37.25,x2=41.75,x3=44.5时,y的数据缺失,将x1、x2、x3的值带入拟合多项式，算出缺失值y1=3.3257,y2=2.0437,y3=4.6002，即可补充缺失数据。4.2 多元线性回归模型1当有缺失的一组数据存在多个自变量时，可以考虑使用多元线性回归模型。我们将data.x

7、ls（见附表一）中的数据5导入matlab(该模型matlab程序见附录二)。首先作出散点图，设定y(PM10)与x1(SO2)、x2(NO2)的关系为二元线性回归模型，即y=b0+b1x1+b2x2。之后作多元回归，求出系数b0=0.03069，b1=-0.023695，b2=0.93619，所以多元线性回归多项式为： Y=0.03069-0.023695*x1+0.93619*x2。且b0、b1、b2在置信度为%95的情况下的置信区间分别为0.024906, 0.036474, -0.08381, 0.036419, 0.76259, 1.1098。再作出残差分析图验证拟合效果，残差较小，

8、说明回归多项式与源数据吻合得较好。若x1=0.055，x2=0.025时，y的数据缺失，则将x1，x2带入回归多项式，算出缺失值y=0.052792。类似地，若x1=0.110,x2=0.045时，y的数据缺失，则将x1，x2带入回归多项式，算出缺失值y=0.070212，即可补充缺失数据。散点图如下：图2 散点图残差分析图如下：图3 残差分析图 五、对问题二的分析和处理5.1 异常值t检验法原理对于问题二，在数据完整但出现异常的情况下，可以考虑使用异常值检验中标准差未知的t检验法。将除可疑测定值以外的其余测定值当做一个总体，并假设该总体服从正态分布。由这些测定值计算平均值与标准差s，而将可疑

9、值当做一个样本容量为1的特殊总体。如果与其余测定值同属于一个总体，则它与其余测定值之间不应有显著性差异。检测统计量为：，假设可由标准差s替代来进行检验,则检测统计量可视为：。若统计量值大于相应置信度下的t检验法的临界值（该临界值通过查表法得出），则将判为异常值2。5.2 举例建模找出异常值我们将data.txt（见附表二）中的数据5导入matlab，首先利用matlab绘制火柴棒图（该模型程序见附录三）观察可疑数据，火柴棒图如下：图4 火柴棒图通过观察上图，得到可疑值为第6,9,13,23,26,29,35,36,45,53行的数据。剔除这些数据后，用剩余数据计算平均值与标准差s。检测统计量为

10、：。置信度=%95，样本容量n=44时，利用查表法易得t检验法的临界值=1.6802。部分t分布表3如下：表2 t分布表n=0.250.100.05410.68051.30251.6829420.68041.30201.6820430.68021.30161.6811440.68011.30111.6802450.68001.30061.6794若统计量值大于临界值，则将判为异常值，根据程序结果，第6,9,13,23,26,29,35,36,45,53行的数据都是异常值。6、对问题三的分析和处理对于问题三，我们采用了分段线性插值、三次样条函数插值以及分段三次Hermite插值法来修正数据异常。

11、同时也需利用外插法修正最后一个数据的异常。本问题使用的是问题二中剔除异常点的数据（数据见附表三）。6.1 原理简介（1）分段线性插值4将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记作，它满足，且在每个小区间上是线性函数。可以表示为有良好的收敛性，即对于有，用 计算x 点的插值时，只用到x 左右的两个节点，计算量与节点个数n 无关。但n 越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。（2）三次样条函数插值4数学上将具有一定光滑性的分段多项式称为样条函数。三次样条函数为：对于上的

12、分划：=b，则，利用样条函数进行插值，即取插值函数为样条函数，称为样条插值。（3）分段三次Hermite插值4如果对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值，这就是Hermite 插值问题。Hermite插值多项式为：6.2 模型求解利用matlab编程（程序见附录四），可作出利用三种插值法修正异常值的图。图如下： 图5 分段线性插值图 图6 三次样条函数插值图 图7 分段三次Hermite插值图根据程序运行结果，可得到分段线性插值的修正结果如表3： 表3 分段线性插值修正结果行号6913232629364553修正结果0.448210.53

13、2290.409490.426830.495380.464510.281250.412640.55435三次样条函数插值的修正结果如表4： 表4 三次样条函数插值修正结果行号6913232629364553修正结果0.426820.594170.432220.432620.519290.432510.433270.410150.28295之后利用第二问中的t检验判断修正后的值是否仍为异常值（matlab程序见附录五）。检验结果显示：分段线性插值和三次样条函数插值所修正的数据中，第36行的数据0.28125和0.43327仍为异常值，第53行的数据0.55435和0.28295也仍为异常值。分段

14、三次Hermite插值的修正结果如表5： 表5 分段三次Hermite插值修正结果行号6913232629364553修正结果0.446890.537940.409490.426830.504190.462640.436180.412640.4689 经过t检验，分段三次Hermite插值所得到的修正值全都不是异常值。所以运用分段三次Hermite插值得到的结果较准确。七、模型评价和推广 问题一采用的一元多项式回归和多元线性回归能够较好地解决补充缺失数据的问题，并且该两种模型兼顾了只有一个自变量和有多个自变量的情况。一元多项式回归均方根误差较小，说明拟合曲线与源数据吻合得较好，故补充的缺失数据

15、准确度也较高。多元线性回归残差较小，说明回归多项式与源数据吻合得较好，故补充的缺失数据是可信的。问题二采用的异常值检验中标准差未知的t检验法能够较准确地判断出异常值。但如果样本容量太小，由标准差s替代来进行检验，可能会产生误差。问题三采用的分段三次Hermite插值能够较好地修正异常数据。并且我们把经过修正的数据采用问题二的方法进行检验，证明所得到的修正值全都不是异常值。因此，修正后的值是较可信的。 缺失数据地补充及异常数据地修正在很多问题中都是必不可少的和基础的，不作这些数据处理就无法解决问题。本文中的模型可以运用到许多需要补充缺失数据或修正异常数据的问题中去。同时，也可利用该模型算法拓展模

16、型在其他方面的适用范围。8、参考文献1 李柏年，MATLAB数据分析方法，北京：机械工业出版社，2012年。2 何正风，MATLAB概率与数理统计分析，北京：机械工业出版社，2012年。3 盛骤，概率论与数理统计，北京：高等教育出版社，2001年。4 司守奎，数学建模算法与程序，北京：国防工业出版社，2007年。5 郑向东，中国大气本底观测元数据，sci-,2013年7月21 日。 9、附录附录一（一元多项式回归模型matlab程序）：clear;x=37.0,37.5,38.0,38.5,39.0,39.5,40.0,40.5,41.0,41.5,42.0,42.5,43.0;y=3.40,

17、3.00,3.00,2.27,2.10,1.83,1.53,1.70,1.80,1.90,2.35,2.54,2.90;plot(x,y,r*)hold onp=polyfit(x,y,2);z=polyval(p,x);plot(x,z,g-);xlabel(x);ylabel(y);legend(实际散点图,拟合曲线);RMSE1=(z-y).2;RMSE=sqrt(sum(RMSE1,2)/13);text(39.0,3.30,拟合曲线为:z=0.16599*x2-13.387*x+271.62);text(39.0,3.50,均方根误差为:RMSE=0.13931);hold off

18、%若x1=37.25,x2=41.75,x3=44.5时,y的数据缺失，则可用得得出的拟合曲线来计算z=0.16599*x2-13.387*x+271.62x1=37.25,x2=41.75,x3=44.5;z1=0.16599*x12-13.387*x1+271.62;z2=0.16599*x22-13.387*x2+271.62;z3=0.16599*x32-13.387*x3+271.62;%计算可得到z1=3.3257,z2=2.0437,z3=4.6002.此例说明数据缺失时可先观察数据，拟合出多项式曲线，再求出缺失的数据。附录二（多元线性回归模型matlab程序）：clear;da

19、ta1=xlsread(H:data.xls);%做出散点图figure(1)scatter3(data1(:,6),data1(:,7),data1(:,8),r.);x=ones(262,1),data1(:,6),data1(:,7);y=data1(:,8);b,bint,r,rint,stats=regress(y,x);xlabel(SO2(X1);ylabel(NO2(X2);zlabel(PM10(Y);text(0.1,0.06,0.2,回归关系式为:Y=0.03069-0.023695*X1+0.93619*X2,color,b);title(PM10（Y）和SO2（X1）

20、、NO2（X2）的关系,color,m);%做残差分析图figure(2)rcoplot(r,rint);xlabel(数据);ylabel(残差);title(残差绘制图);%补缺失数据x1=1 0.055 0.025;y1=x1*b;x2=1 0.110 0.045;y2=x2*b;%y1=0.052792;y2=0.070212. 附录三（t检验法matlab程序）：clear;load H:data.txt;load H:data1.txt;s=sum(data(:,7),1)/53;alpha=0.05;n=44;T=1.6802;y=1;a=1;stem(data(:,7);%观察

21、异常值.%观察异常值.可得到异常值为第6,9,13,23,26,29,35,36,45,53行.%把不是异常的个体看成一个整体data1,异常个体看成data2data2=1.004571;0.888936;0.77689;0.856265;0.79843;0.935673;0.630218;0.812906;0.101234;aver=sum(data1(:,7),1)/44;s=sqrt(sum(data1(:,7)-aver).2,1)/44);for y=1:9 K(y,1)=abs(data2(y,1)-aver)./s; if(K(y,1)T) H(a,1)=y; a=a+1; e

22、nd end 附录四（三种插值法修正异常值程序）：clear;load H:data1.txt;x0=1,2,3,4,5,7,8,10,11,12,14,15,16,17,18,19,20,21,22,24,25,. 27,28,30,31,32,33,34,35,37,38,39,40,41,42,43,44,46,47,48,49,50,51,52;y0= 0.442955 0.39876 0.37345 0.437790 0.44321 0.45321 0.54312 0.52145 0.34582 0.343307. 0.475664 0.423589 0.432951 0.44124

23、6 0.421128 0.416297 0.435372 0.394432 0.428593. 0.425061 0.461622 0.529131 0.460915 0.468112 0.523575 0.490494 0.472183 0.483644. 0.456781 0.425717 0.431464 0.432782 0.385745 0.388722 0.401138 0.385566 0.423361 0.401917. 0.443387 0.445423 0.573121 0.436781 0.447899 0.501123; x=0:0.5:52; y1=interp1(x

24、0,y0,x); y2=interp1(x0,y0,x,spline);pp1=csape(x0,y0); y3=ppval(pp1,x);y6=interp1(x0,y0,x,pchip);y5=interp1(x0,y0,53,spline,extrap);y4=interp1(x0,y0,53,linear,extrap);y7=interp1(x0,y0,53,pchip,extrap);subplot(2,2,1), plot(x0,y0,+,x,y1), title(Piecewise linear); subplot(2,2,2), plot(x0,y0,+,x,y2), tit

25、le(Spline);subplot(2,2,3), plot(x0,y0,+,x,y6), title(pchip);y1=y1;y2=y2;y3=y3;y6=y6; % x=53时，y4=0.55435,y5=0.28292,y7=0.4689附录五（修正后的值作t检验的程序）：clear;load H:data1.txt;alpha=0.05;n=44;T=1.6802;y=1;a=1;%观察异常值.可得到异常值为第6,9,13,23,26,29,36,45,53行.%把不是异常的个体看成一个整体data1,异常个体看成data3%用线性插值data3=0.44821;0.53229;0

26、.40949;0.42683;0.49538;0.46451;0.28125;0.41264;0.55435;aver=sum(data1(:,7),1)/44;s=sqrt(sum(data1(:,7)-aver).2,1)/44);for y=1:9 K1(y,1)=abs(data3(y,1)-aver)./s; if(K1(y,1)T) H1(a,1)=y; a=a+1; end end%用三次样条插值b=1;data4=0.42682;0.59417;0.43222;0.43262;0.51929;0.43251;0.43327;0.41015;0.28295; for y=1:9

27、K2(y,1)=abs(data4(y,1)-aver)./s; if(K2(y,1)T) H2(b,1)=y; b=b+1; end end %用分段三次Hermite插值 c=1; data5=0.44689;0.53794;0.40949;0.42683;0.50419;0.46264;0.43618;0.41264;0.4689; for y=1:9 K3(y,1)=abs(data5(y,1)-aver)./s; if(K3(y,1)T) H3(c,1)=y; c=c+1; end end 附表一（多元线性回归模型使用的数据）：城市名称年月日点位代码点位名称SO2NO2PM10A城2

28、0041210110.1160.0740.074A城20041310110.1200.0780.075A城20041410110.0620.0530.048A城20041510110.0500.0420.021A城20041610110.1310.0610.054A城20041710110.0710.0350.025A城20041810110.0480.0350.012A城20041910110.0450.0350.011A城200411110110.1040.0500.032A城200411210110.1500.0500.037A城200411410110.0750.0510.065A城200411510110.0250.0280.032A城200411610110.0170.0280.037A城200411710110.0170.0260.052A城200411810110.0300.0300.077A城200412010110.0570.0300.076A城200412110110.0530.0170.062A城200412210110.0200.0230.064A城2004