缺失数据地补充及异常数据地修正Word格式.docx
《缺失数据地补充及异常数据地修正Word格式.docx》由会员分享，可在线阅读，更多相关《缺失数据地补充及异常数据地修正Word格式.docx（39页珍藏版）》请在冰点文库上搜索。
来进行检验。
3、符号说明
符号
含义
RMSE
均方根误差
样本平均值
s
样本标准差
可疑测定值
总体标准差
k
检测统计量
置信度
下的t检验法临界值
n
样本容量
4、对问题一的分析和处理
4.1一元多项式回归模型[1]
当有缺失的一组数据只有一个自变量时，可以考虑使用一元多项式回归模型。
我们采用的数据如下表：
表1
首先使用matlab作出散点图，观察其形状，决定拟合多项式的次数，运用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.xls（见附表一）中的数据[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]。
再作出残差分析图验证拟合效果，残差较小，说明回归多项式与源数据吻合得较好。
散点图如下：
图2散点图
残差分析图如下：
图3残差分析图
五、对问题二的分析和处理
5.1异常值t检验法原理
对于问题二，在数据完整但出现异常的情况下，可以考虑使用异常值检验中标准差未知的t检验法。
将除可疑测定值
与标准差s，而将可疑值
当做一个样本容量为1的特殊总体。
如果
与其余测定值同属于一个总体，则它与其余测定值之间不应有显著性差异。
（该临界值通过查表法得出），则将
判为异常值[2]。
5.2举例建模找出异常值
我们将data.txt（见附表二）中的数据[5]导入matlab，首先利用matlab绘制火柴棒图（该模型程序见附录三）观察可疑数据，火柴棒图如下：
图4火柴棒图
通过观察上图，得到可疑值为第6,9,13,23,26,29,35,36,45,53行的数据。
剔除这些数据后，用剩余数据计算平均值
与标准差s。
=%95，样本容量n=44时，利用查表法易得t检验法的临界值
=1.6802。
部分t分布表[3]如下：
表2t分布表
α=0.25
0.10
0.05
41
0.6805
1.3025
1.6829
42
0.6804
1.3020
1.6820
43
0.6802
1.3016
1.6811
44
0.6801
1.3011
1.6802
45
0.6800
1.3006
1.6794
若统计量值大于临界值
，则将
判为异常值，根据程序结果，第6,9,13,23,26,29,35,36,45,53行的数据都是异常值。
6、对问题三的分析和处理
对于问题三，我们采用了分段线性插值、三次样条函数插值以及分段三次Hermite插值法来修正数据异常。
本问题使用的是问题二中剔除异常点的数据（数据见附表三）。
6.1原理简介
（1）分段线性插值[4]
将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记作
，它满足
，且
在每个小区间
上是线性函数
可以表示为
有良好的收敛性，即对于
有，
用
计算x点的插值时，只用到x左右的两个节点，计算量与节点个数n无关。
但n越大，分段越多，插值误差越小。
实际上用函数表作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。
（2）三次样条函数插值[4]
数学上将具有一定光滑性的分段多项式称为样条函数。
三次样条函数为：
对于
上的分划
：
=b，则，
利用样条函数进行插值，即取插值函数为样条函数，称为样条插值。
（3）分段三次Hermite插值[4]
如果对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值，这就是Hermite插值问题。
Hermite插值多项式为：
6.2模型求解
利用matlab编程（程序见附录四），可作出利用三种插值法修正异常值的图。
图如下：
图5分段线性插值图图6三次样条函数插值图
图7分段三次Hermite插值图
根据程序运行结果，可得到分段线性插值的修正结果如表3：
表3分段线性插值修正结果
行号
6
9
13
23
26
29
36
53
修正
结果
0.44821
0.53229
0.40949
0.42683
0.49538
0.46451
0.28125
0.41264
0.55435
三次样条函数插值的修正结果如表4：
表4三次样条函数插值修正结果
0.42682
0.59417
0.43222
0.43262
0.51929
0.43251
0.43327
0.41015
0.28295
之后利用第二问中的t检验判断修正后的值是否仍为异常值（matlab程序见附录五）。
分段线性插值和三次样条函数插值所修正的数据中，第36行的数据0.28125和0.43327仍为异常值，第53行的数据0.55435和0.28295也仍为异常值。
分段三次Hermite插值的修正结果如表5：
表5分段三次Hermite插值修正结果
修正结果
0.44689
0.53794
0.50419
0.46264
0.43618
0.4689
经过t检验，分段三次Hermite插值所得到的修正值全都不是异常值。
七、模型评价和推广
问题一采用的一元多项式回归和多元线性回归能够较好地解决补充缺失数据的问题，并且该两种模型兼顾了只有一个自变量和有多个自变量的情况。
一元多项式回归均方根误差较小，说明拟合曲线与源数据吻合得较好，故补充的缺失数据准确度也较高。
多元线性回归残差较小，说明回归多项式与源数据吻合得较好，故补充的缺失数据是可信的。
问题二采用的异常值检验中标准差未知的t检验法能够较准确地判断出异常值。
但如果样本容量太小，由标准差s替代σ来进行检验，可能会产生误差。
问题三采用的分段三次Hermite插值能够较好地修正异常数据。
并且我们把经过修正的数据采用问题二的方法进行检验，证明所得到的修正值全都不是异常值。
因此，修正后的值是较可信的。
缺失数据地补充及异常数据地修正在很多问题中都是必不可少的和基础的，不作这些数据处理就无法解决问题。
本文中的模型可以运用到许多需要补充缺失数据或修正异常数据的问题中去。
同时，也可利用该模型算法拓展模型在其他方面的适用范围。
8、参考文献
[1]李柏年，MATLAB数据分析方法，北京：
机械工业出版社，2012年。
[2]何正风，MATLAB概率与数理统计分析，北京：
[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,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*'
）
holdon
p=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*x^2-13.387*x+271.62'
text（39.0,3.50,'
均方根误差为:
RMSE=0.13931'
holdoff
%若x1=37.25,x2=41.75,x3=44.5时,y的数据缺失，则可用得得出的拟合曲线来计算z=0.16599*x^2-13.387*x+271.62
x1=37.25,x2=41.75,x3=44.5;
z1=0.16599*x1^2-13.387*x1+271.62;
z2=0.16599*x2^2-13.387*x2+271.62;
z3=0.16599*x3^2-13.387*x3+271.62;
%计算可得到z1=3.3257,z2=2.0437,z3=4.6002.此例说明数据缺失时可先观察数据，拟合出多项式曲线，再求出缺失的数据。
附录二（多元线性回归模型matlab程序）：
data1=xlsread（'
H:
\data.xls'
%做出散点图
figure
（1）
scatter3（data1（:
6）,data1（:
7）,data1（:
8）,'
r.'
x=[ones（262,1）,data1（:
7）];
y=data1（:
8）;
[b,bint,r,rint,stats]=regress（y,x）;
SO2（X1）'
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）、NO2（X2）的关系'
m'
%做残差分析图
figure
（2）
rcoplot（r,rint）;
数据'
残差'
残差绘制图'
%补缺失数据
x1=[10.0550.025];
y1=x1*b;
x2=[10.1100.045];
y2=x2*b;
%y1=0.052792;
y2=0.070212.
附录三（t检验法matlab程序）：
loadH:
\data.txt;
\data1.txt;
s=sum（data（:
7）,1）/53;
alpha=0.05;
n=44;
T=1.6802;
y=1;
a=1;
stem（data（:
7））;
%观察异常值.
%观察异常值.可得到异常值为第6,9,13,23,26,29,35,36,45,53行.
%把不是异常的个体看成一个整体data1,异常个体看成data2
data2=[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）;
fory=1:
K（y,1）=abs（data2（y,1）-aver）./s;
if（K（y,1）>
T）
H（a,1）=y;
a=a+1;
end
end
附录四（三种插值法修正异常值程序）：
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.4429550.398760.373450.4377900.443210.453210.543120.521450.345820.343307...
0.4756640.4235890.4329510.4412460.4211280.4162970.4353720.3944320.428593...
0.4250610.4616220.5291310.4609150.4681120.5235750.4904940.4721830.483644...
0.4567810.4257170.4314640.4327820.3857450.3887220.4011380.3855660.4233610.401917...
0.4433870.4454230.5731210.4367810.4478990.501123];
x=0:
0.5:
52;
y1=interp1（x0,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,'
extrap'
y4=interp1（x0,y0,53,'
linear'
y7=interp1（x0,y0,53,'
subplot（2,2,1）,plot（x0,y0,'
+'
x,y1）,title（'
Piecewiselinear'
subplot（2,2,2）,plot（x0,y0,'
x,y2）,title（'
Spline'
subplot（2,2,3）,plot（x0,y0,'
x,y6）,title（'
y1=y1'
;
y2=y2'
y3=y3'
y6=y6'
%x=53时，y4=0.55435,y5=0.28292,y7=0.4689
附录五（修正后的值作t检验的程序）：
%观察异常值.可得到异常值为第6,9,13,23,26,29,36,45,53行.
%把不是异常的个体看成一个整体data1,异常个体看成data3
%用线性插值
data3=[0.44821;
0.53229;
0.40949;
0.42683;
0.49538;
0.46451;
0.28125;
0.41264;
0.55435];
K1（y,1）=abs（data3（y,1）-aver）./s;
if（K1（y,1）>
H1（a,1）=y;
%用三次样条插值
b=1;
data4=[0.42682;
0.59417;
0.43222;
0.43262;
0.51929;
0.43251;
0.43327;
0.41015;
0.28295];
fory=1:
K2（y,1）=abs（data4（y,1）-aver）./s;
if（K2（y,1）>
H2（b,1）=y;
b=b+1;
%用分段三次Hermite插值
c=1;
data5=[0.44689;
0.53794;
0.50419;
0.46264;
0.43618;
0.4689];
K3（y,1）=abs（data5（y,1）-aver）./s;
if（K3（y,1）>
H3（c,1）=y;
c=c+1;
附表一（多元线性回归模型使用的数据）：
城市名称
年
月
日
点位代码
点位名称
SO2
NO2
PM10
A城
2004
1
2
101
0.116
0.074
3
0.120
0.078
0.075
4
0.062
0.053
0.048
5
0.050
0.042
0.021
0.131
0.061
0.054
7
0.071
0.035
0.025
8
0.012
0.045
0.011
11
0.104
0.032
12
0.150
0.037
14
0.051
0.065
15
0.028
16
0.017
17
0.026
0.052
18
0.030
0.077
20
0.057
0.076
21
22
0.020
0.023
0.064
升级会员 