题目来源于2019 全国大学生数学建模大赛D题-空气质量数据校准。 通过对两尘四气的危害，空气污染对生态环境和人类健康危害极大（PM2.5、PM10、CO、NO2、SO2、O3)实时监测浓度可及时掌握空气质量，对污染源采取相应措施。虽然国家监测控制站(国家控制点)对两尘四气有监测数据，比较准确，但由于国家控制点的控制较少，数据发布时间滞后较长，成本较高，无法对实时空气质量进行监测和预测。公司自主开发的微型空气质量探测器（如图所示）成本低，可实时监测某一地区的空气质量，监测温度、湿度、风速、气压、降水等气象参数。

本次实战数据分析涉及以下技术内容： (1)处理重复值 (2)时间类型数据转换 （3）插值方法 （4）归一化 （5）标准化 (6)灰色关联 (7)回归模型的培训和评价 (8)评价各自变量的重要性 （9）NCA算法选择特征 (10)新数据集预测模型

以下主要用途MATLAB功能实现： （1）readtable （2）unique （3）ismissing （4）table2array （5）datenum （6）interp1 （7）mapminmax （8）mapstd （9）fitrauto （10）fsrnca （11）setdiff （12）predict

读取数据，查看数据的基本情况，

warning("off"); Data_1=readtable("附件1.csv"); Data_2=readtable("附件2.csv"); size(Data_1) size(Data_2) 

输出

由于后续的插值和回归都是函数问题，因此需要确保一个时间戳只能对应一个y值。

[u_1,ia_1,ic_1]=unique(Data_1{:,end}); [u_2,ia_2,ic_2]=unique(Data_2{:,end}); u_Data_1=Data_1(ia_1,:); u_Data_2=Data_2(ia_2,:); size(u_Data_1) size(u_Data_2) 

输出

使用函数readtable之后，缺失值NaN形式存在，

miss_value_1=sum(sum(ismissing(u_Data_1)==true)) miss_value_2=sum(sum(ismissing(u_Data_2)==true)) 

输出 国建点和自建点数据均无缺失值。

因为我不是该项目涉及的专业领域的专家，我没有定义异常值的经验，我认为统计意义上的异常并不意味着特定问题中的数据确实异常，考虑到没有0值得存在，所以我没有在这里处理数据。

%使用国建点时间戳插入自建点数据 time1=u_Data_1{:,end}; time2=u_Data_2{:,end}; time1_num=datenum(time1)；%时间数据转换为数字型 time2_num=datenum(time2)；将%时间类型的数据转换为数字类型  for i=1:6 %"两尘四气"共6列     Data_1_spline(:,i)=interp1(time2_num,u_Data_2{:,i},time1_num,"spline");     Data_1_linear(:,i) = interp1(time2_num,u_Data_2{:,i},time1_num,"linear"); end 

输出

对比PM2.5.对于插值前后的性能，样条插值对时间间隔大、变化陡峭点的插值敏感，容易出现异常插值，线性插值性能相对较好。

虽然线性插值表现良好，但为了进一步提高插值的准确性，我计划对每个相邻时间戳的间隔进行统计分析，以设定时间间隔的边界值。如果间隔大于此值，我会认为插值误差大，可信度低；相反，我认为可信度高，插值可取。 分析自建点的时间间隔，

%分析自建点数据中两个时间点的间隔 jiange_time2_num=[]; for i=2:length(time2_num)     jiange_time2_num(i-1)=time2_num(i)-time2_num(i-1); end  jiange_time2_num_t=tabulate(jiange_time2_num); jiange=jiange_time2_num_t(:,1); bei=[]; for i=2:length(jiange)     bei(i-1)=jiange(i)/jiange(1)min倍数，也就是分钟 end 

时间间隔频率表如下，间隔最小，频率最大为1分钟。

以上，为了保证后续插值的有效性，我只在一定范围内插入时间间隔小于国家建设点的时间戳在这里确定间隔是30min为可接受的有效插值间隔。计算此时有效插值间隔的比值，

sum_t_57=sum(jiange_time2_num_t(1:57,3)) 

得到 发现有效插值间隔占总时间间隔99.9492%几乎包含了大部分时间点。因此，我计划逐一寻找每一个time1_num(k)在time2_num里的位置m~m 1，并计算time2_num(m)到time2_num(m 1)距离，如果距离小于30min，则留下该点time1_num(k)进行插值。

%确定time1_num有效插值的点time1_num的 %第一点在time2_num在第一点之前，所以从2开始 liu=[]; for i=2:length(time1_num)     ind_1=find(time2_num<time1_num(i));     ind_2=ind_1(end);     if (time2_num(ind_2 1)-time2_num(ind_2)) <0.0208         liu(i)=1;     end   end liu(1)=1;  sum(liu==1)/length(liu) 

筛选后，计算得出time1_num保留的时间戳数量占time1_num总数的96.3%，绝大多数点被选中。

筛选后使用time1_num重新插入两尘四气等五个影响因素，并绘图比较。spline方法容易出现异常点，这里只用linear方法。

time1_num_new=time1_num(liu==1); Data_1_spline=[]; Data_1_linear=[]; for i=1:6 %"两尘四气"共6列     Data_1_linear(:,i) = interp1(time2_num,u_Data_2{:,i},time1_num_new,"linear"); end  j=1; Data_1linear_factor=[];

for i=7:11 %影响因素共5列
    Data_1_linear_factor(:,j) = interp1(time2_num,u_Data_2{:,i},time1_num_new,"linear");
    j=j+1;
end
Data_1_linear_factor(1,:) =Data_1_linear_factor(2,:);
figure,
plot(time2_num,u_Data_2{:,4},'-ro',time1_num_new,Data_1_linear(:,4),'bo');
figure,plot(time2_num,u_Data_2{:,4},'-ro');
figure,plot(time1_num_new,Data_1_linear(:,4),'bo');

对比图， 可见由于我们设定了最小插值时间间隔，插值函数未对间隔大于30min的时间戳进行插值，尽量保证了插值的准确度，避免因为插值产生的误差在后续求关联度和进行回归时被放大。

先进行灰色关联前数据预处理。网上关于这个预处理，大家用的方法各有不同。有归到同一起点的，有归到平均值的（除以平均值），有归到特定范围内的（如0~1），而且算出来的关联度都互不相同。我查阅网上大家对灰色关联的解释，觉得既然灰色关联是表征两组曲线变化趋势相似度的，那归到固定范围，应该最为合适。因为我的理解，趋势即为斜率，斜率即为比例，既然只考虑比例，那大家都归到同一固定范围不就好了？所以以下我采用归到固定范围内来预处理各列数据。如果有数学大神，请给予我理论性的指导，谢谢。

Data_1_num=table2array(u_Data_1(liu==1,1:end-1));
error_data=Data_1_linear-Data_1_num;
%误差归一化
error_data_B=mapminmax(error_data',0,1);
error_data_B=error_data_B';
error_data_B(1,:)=error_data_B(2,:);
%自变量归一化
Data_X=Data_1_linear_factor;
Data_X_B=mapminmax(Data_X',0,1);
Data_X_B=Data_X_B';

灰色关联算法求关联度，

GLD=[];
for i=1:size(error_data_B,2)
    GLD(i,:)=GRA2(Data_X_B,error_data_B(:,i));  
end
GLD

关联度结果表， 所使用的灰色关联函数如下，

function out=GRA2(X,Y)
absX0_Xi=abs(X-repmat(Y,1,size(X,2)));
a=min(min(absX0_Xi));
b=max(max(absX0_Xi));
rho=0.5;
gamma=(a+rho*b)./(absX0_Xi+rho*b);
out=mean(gamma);
end

先进行所需数据预处理。首先是数据标准化，当样本量足够大的时候，随机划分的训练集和测试集其实是有相同的分布的，这里在划分数据集前进行数据标准化，

[y_data,ps_1]=mapstd(error_data');%因为要向训练集范围外回归，所以未使用归一化
y_data=y_data';

[x_data,ps_2]=mapstd(Data_X');
x_data=x_data';

构建回归模型所需数据，先取误差的第一列，即将自建点和国建点PM2.5的差值作为回归目标，进行训练，

tbl_1=array2table([y_data(:,1),x_data]);
tbl_1.Properties.VariableNames=["Y","x1","x2","x3","x4","x5"];

划分数据集，

rng('default') % For reproducibility of the data partition
c = cvpartition(size(tbl_1,1),'Holdout',0.2);
trainingIdx = training(c); % Training set indices
tbl_1_Train = tbl_1(trainingIdx,:);
testIdx = test(c); % Test set indices
tbl_1_Test = tbl_1(testIdx,:);

训练回归模型，使用了fitrauto，该函数可以自动选择最优模型及对应的超参数，是我等懒人“居家旅行，摸鱼放松”的必备函数，(•‾̑⌣‾̑•)

options = struct('UseParallel',true);
Mdl = fitrauto(tbl_1_Train,'Y','OptimizeHyperparameters','all', ...
    'HyperparameterOptimizationOptions',options,...
    "Learners","tree");

这里我只选定使用树模型，便于后续对各特征打分。大家也可以尝试不同的模型，或者采用大人的方法：“我都要！”

评价模型在测试集上的性能，

testMSE = loss(Mdl,tbl_1_Test,'Y');
testError = log(1 + testMSE)

模型对各自变量在回归模型中的重要性进行打分，

predictorImportance(Mdl)

使用NCA算法做特征选择，也可以得出各自变量在回归时的重要性排序，

X=tbl_1_Train{:,2:end};
y=tbl_1_Train{:,1};
N=length(y);
mdl_nca = fsrnca(X,y,'Verbose',1,'Lambda',0.5/N);
mdl_nca.FeatureWeights

可以看出，两种对特征重要度的评价方法都对第一个特征打了最低分。我觉得可以结合上述灰色关联度，对各影响因素进行综合评价。

使用训练好的模型在新数据上做预测，这里相当于“外推”，向外拟合，

[~,ia] = setdiff(time2_num,time1_num_new);
y_base=u_Data_2{ia,1};
predict_x=u_Data_2{ia,7:11};
predict_x=mapstd('apply',predict_x',ps_2);
predict_x=predict_x';
Yfit = predict(Mdl,predict_x);
%逆标准化
Yfit = mapstd('reverse',Yfit',ps_1);
Yfit =Yfit';
y_new=y_base-Yfit;

画图比对，

figure,
plot(time1_num,u_Data_1{:,1},'ro','LineWidth',2),hold on;
plot(time2_num(ia),y_new(:,1),'bo');hold off;

figure,
plot(time1_num,u_Data_1{:,1},'ro','LineWidth',2),hold on;
plot(time2_num,u_Data_2{:,1},'go');hold off;

可见，未修正前的自建点数据（图二），整体要比国建点偏上，而修正后的自建点数据（图一），整体将国建点“包裹”的较好，更加合理。当然，使用回归的方法修正后，自建点出现了一些负值，将这些负值置0即可。 大家可以尝试对剩下的“一尘四气”进行回归，并比对效果，此处不赘述了。