MATLAB曲⾯插值及交叉验证
在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。插值是离散函数逼近的重要⽅法，利⽤它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。曲⾯插值是对三维数据进⾏离散逼近的⽅法，MATLAB中的曲⾯插值函数有Triscatteredinterp，interp2，griddata等。我们以griddata为例讲解曲⾯插值及其交叉验证的过程。
⼀、  gridata曲⾯插值
gridata不仅可以对三维曲⾯进⾏插值，还能对四维的超平⾯进⾏插值。griddata的调⽤⽅法：ZI = griddata(x,y,z,XI,YI,method);[XI,YI,ZI] =
griddata(x,y,z,XI,YI,method)。插值⽅法有'linear'以三⾓形为基础的线性内插、'cubic'以三⾓形为基础的三次⽅程内插，'natural'⾃然邻近插值，'nearest'最近邻近插值，'v4'由MATLAB提供的插值⽅法，默认值为'linear'。三⾓形为基础，就是按Delaunay⽅法先出内插点四周的3个点，构成三⾓形，内插点在三⾓形内。然后线性内插，或三次⽅程内插。 'cubic' 和 'v4' 插值结果构成的曲⾯较光滑。
我们可以⽤griddata函数来求得未知点的Z值，进⽽绘制曲⾯的⼤致图像，代码如下
[X,Y]=meshgrid(min(x):0.1:max(x),min(y):0.1:max(y)); %确定⽹络坐标
Z=griddata(x,y,z,X,Y,'cubic');              %在⽹格点位置插值求Z
surf(X,Y,Z)                                  %绘制曲⾯
title 'Points to Surface by griddata'
所得的结果如下
⼆、拟合精度指标
做完曲⾯插值后，如何对曲⾯插值模型的拟合好坏进⾏评判是我们应该考虑的问题。MATLAB曲⾯拟合⼯具箱cftool中⽆论曲⾯插值还是函数拟合都⽤拟合优度R-squared来评判拟合效果的好坏。其公式为
其中，y为实际值， y^为预测值， y*为y的平均值。
个⼈觉得函数拟合可以⽤此指标评判，但曲⾯插值⽤此指标有些不妥。故介绍⼀个新的评判插值效果的指标，公式为
符号含义与上式相同，虽然R-squared有很多很好的数学性质，但对于前提假设不明确的模型来说，⽤拟合度指标r显得更直观些，即1扣除残差平均占实际平均的⽐例。
下⾯我们来对⽐两个指标随样本点增加的变化情况，代码如下
function r=FitLevel(y0,y1) %拟合度指标的公式。y0为真实值，y1为预测值
%r=sum((y1-mean(y0)).^2)/sum((y0-mean(y0)).^2);%拟合优度R-squared
r=1-sqrt(sum((y1-y0).^2)/sum(y0.^2)); %拟合度指标
%绘制拟合度指标随样本数增加的变化图
step=1;
I=1000:1000:length(x)
for i=I
Z=griddata(x(1:i),y(1:i),z(1:i),x(1:i),y(1:i));
r=FitLevel(z(1:i),Z); %求拟合度指标
R(step)=r;
step=step+1
end
plot(I,R)
title '拟合度-样本点数'
拟合优度和新的拟合度随样本点增加的变化情况图如下
虽说这两幅图曲线的形状有点类似，但拟合优度那副⽐较平稳，新的拟合度波动幅度较⼤，从数值上分析，还是新的拟合度指标较为合理的。在交叉验证那节还可发现⽤拟合优度去检验预测效果值⼤于1，超出它的取值范围，故若要⽤交叉验证算法检验插值效果还是得⽤新的拟合度指标的。
三、交叉验证
常⽤的精度测试⽅法主要是交叉验证，交叉验证的基本思想是把在某种意义下将原始数据(dataset)进⾏分组,⼀部分做为训练集(train set)，另⼀部分做为验证集
(validation set or test set),⾸先⽤训练集对分类器进⾏训练，再利⽤验证集来测试训练得到的模型(mode
l)，以此来做为评价分类器的性能指标。在此基础上，反复的做训练、测试以及模型的选择。包含简单交叉验证，K折交叉验证，留⼀验证，Holdout 验证等。结合业务场景，本⽂将把K折交叉验证和简单交叉验证结合使⽤。      K折交叉验证，初始采样分割成K个⼦样本，⼀个单独的⼦样本被保留作为验证模型的数据，其他K-1个样本⽤来训练。交叉验证重复K次，每个⼦样本验证⼀次，平均K次的结果或者使⽤其它结合⽅式，最终得到⼀个单⼀估测。这个⽅法的优势在于，同时重复运⽤随机产⽣的⼦样本进⾏训练和验证，每次的结果验证⼀次。10折交叉验证是最常⽤的。由于本次业务场景的数据并⾮空间数据，⽽是时序数据，所以并没有重复验证K次，⽽是⽤前K-1份样本⽤于训练模型，第K份样本⽤于测试模型。
四、应⽤实例
本次使⽤交叉验证的⽬的主要是为了在已测得⾜够量数据的情况下，寻出适⽤于该样本的最优插值模型和该模型对应的最优样本量。
设训练模型的拟合度为r_train，理想阈值为train；测试模型的拟合度为r_test，理想阈值为test；整体拟合效果为:
r_total=weight*r_train+(1-weight)*r_train
其中weight为训练数据拟合度占整体拟合度重要程度的百分⽐，默认值0.5。求某⼀插值⽅法对应最优或理想样本量的算法如下：
1）建模次数step=1，步长为stepsize。当前训练样本量i
2）对（1,i）的所有数据进⾏插值，求出模型的拟合度r_train。
3）对（i,(K/K-1)*i）的所有数据进⾏检验，求出检测的精度r_test。
4）判断r_total(step)> r_total(step-1)，若是，则记录样本数i。
5）判断r_train>train且r_test>test，若是，则记录样本数i。
6）Step=step+1;i=i+stepsize；返回2)，直到循环结束。
通过以上算法即可求出整体情况最优时的样本数和⼈⼯调节阈值时的样本数序列。根据此算法设计了CrossValidation函数，详细内容参见函数源⽂件，其输⼊输出的完整参数如下：
[Result,R_total]=CrossValidation(x,y,z,train,test,weight,method,stepsize,folder)
我们先输⼊⼀些参数看看函数的输出结果
Result=CrossValidation(x,y,z,0.89,0.7,0.6,'cubic',1000)
最终打印出了训练模型拟合度指标r_total随样本数增加⽽变化的图形、测试数据精度指标r_test随样本数增加⽽变化的图形和总体拟合效果随样本数增加⽽变化的图形。
还返回了⼀个Result矩阵，第⼀⾏为样本数i，第⼆⾏为训练模型拟合度r_train，第三⾏为测试拟合度r_test。第⼀列为总体拟合效果最⼤时返回的值，第⼆列及之后的列是r_train>train且r_test>test时返回的值。可见，275000左右的样本数就能很好的⽤于该插值模型了。
该函数还可以输出整体拟合指标序列R_total，可⽤于对⽐不同⽅法对该样本的拟合⽔平。通过判断sum(R_total1>R_total2)/length(R_total1)是否⼤于0.5来判断谁优谁劣，若⼤于0.5，则说明输出R_total1序列的⽅法拟合⽔平较⾼。下⾯是⼀段对⽐的实例。
[R1,r1]=CrossValidation(x,y,z,0.89,0.7,0.6,'nearest',10000);
[R2,r2]=CrossValidation(x,y,z,0.89,0.7,0.6,'linear',10000);
[R3,r3]=CrossValidation(x,y,z,0.89,0.7,0.6,'natural',10000);
[R4,r4]=CrossValidation(x,y,z,0.89,0.7,0.6,'cubic',10000);
sum(r1>r2)/length(r1)  %返回0，说明r2所有值都⼤于r1
sum(r3>r2)/length(r2)  %返回1，说明r3所有值都⼤于r2
sum(r4>r3)/length(r3)  %返回0，说明r3所有值都⼤于r4
可见，⾃然邻近插值法拟合精度⾼于其他⼏种⽅法，⽐较适⽤于该样本。不过对⽐R1、R2、R3、R4后发现这四种⽅法都能出相同的最优样本点，应该是拟合度指标变化趋势相同的原因。
以上论述已实现求理想样本点和理想插值⽅法的需求。由于我使⽤的⼯程数据具有保密性，不便附上源数据，可⾃⼰随机⽣成数据进⾏测试，下⾯附上对插值模型进⾏交叉验证的函数源码。
function [Result,R_total]=CrossValidation(x,y,z,train,test,weight,method,stepsize,folder)
% CrossValidation：曲⾯插值交叉验证。把样本分为训练数据和测试数据，并对训练数据
%                进⾏曲⾯插值，求出起训练集拟合精度，⽤该模型对测试数据进⾏预
%                测，求出测试集的拟合精度。最终求出或⼈为判断出该样本较为合理
%                的⽤于拟合插值模型的样本数。
%
% usage #1: [Result,R_total]=CrossValidation(x,y,z);
% usage #2: [Result,R_total]=CrossValidation(x,y,z,train,test);
% usage #3: [Result,R_total]=CrossValidation(x,y,z,train,test,weight);
% usage #4: [Result,R_total]=CrossValidation(x,y,z,train,test,weight,method);
% usage #5: [Result,R_total]=CrossValidation(x,y,z,train,test,weight,method,stepsize);
% usage #6: [Result,R_total]=CrossValidation(x,y,z,train,test,weight,method,stepsize,folder);
%
% Arguments: (input)
%x,y,z - 输⼊向量
%train - 训练数据拟合度的理想阈值，范围（0,1），默认值0.8
%test - 测试数据拟合度的理想阈值，范围（0.1），默认值0.7
%weight - 训练拟合度相对于测试拟合度的权重，范围（0.1），默认值0.5
%method - 插值所⽤的⽅法，默认值'lineat'，可选范围如下：
%      'nearest'  - Nearest neighbor interpolation
%      'linear'    - Linear interpolation (default)
%      'natural'  - Natural neighbor interpolation
%      'cubic'    - Cubic interpolation (2D only)
%      'v4'        - MATLAB 4 griddata method (2D only)
%stepsize - 步长，⼦样本每次增加的数⽬，默认值为5000。matlab拟合数据
%folder - K折交叉验证的折数，默认值10.意味着把样本数分为10份，拿1份去验证。
%
%Result：返回⼀个Result矩阵和⼀个整体拟合效果队列R_total；⼀副训练数据随点数增
%        加的拟合效果图R_train-number，⼀副测试数据随点数增加的拟合效果图
%        R_test-number，⼀副整体情况随点数增加的拟合效果图R_total-number。
%
%        Result = [  i(m)              i(n)          ...
%                    R_train(i(m))      R_train(i(n))  ...
%                    R_test(i(m))      R_train(i(n))  ...  ]
%        第⼀列为整体拟合效果最⼤时对应的点数、训练数据拟合度，测试数据拟合度。
%        第⼆列及之后所有列为训练数据⼤于理想训练数据阈值train且测试数据⼤于测
%        试数据理想阈值test时的样本点数、训练数据拟合度，测试数据拟合度。
%
%
%Author: ⼩江同学
%e-mail: 842419386@qq
%设置参数默认值
if nargin==3
train=0.8;test=0.7;weight=0.5;method='linear';folder=10;stepsize=5000;
elseif nargin==5
weight=0.5;method='linear';folder=10;stepsize=5000;
elseif nargin==6
method='linear';folder=10;stepsize=5000;
elseif nargin==7
folder=10;stepsize=5000;
elseif nargin==8
folder=10;
end
%去除输⼊向量的空值
k = isnan(x) | isnan(y) | isnan(z);
if any(k)
x(k)=[];
y(k)=[];
z(k)=[];
end
%变量初始化
step=1;      %循环次数
r_col=2;      %返回Result矩阵的列数
r_max=0;      %r_total的最⼤值
R_train=[];  %训练拟合度序列
R_test=[];    %测试拟合度序列
Result=[];    %返回结果
h1=plot(0,0); %动画初始化
hold on
h2=plot(0,0);
total=weight*train+(1-weight)*test;  %整体拟合效果阈值
I=round(0.3*length(x)):stepsize:round((1-1/folder)*length(x));  %循环范围及步长
%交叉验证及寻最优或理想样本数
for i=I
Z=griddata(x(1:i),y(1:i),z(1:i),x,y,method);  %x(1:i),y(1:i),z(1:i)建模，预测x,y所对应的z值
z_temp=z;
N=find(isnan(Z));            %去除预测向量中的空值
Z(N)=[];z_temp(N)=[];
r_train=FitLevel(z_temp(1:i),Z(1:i));        %求训练模型拟合度
r_test=FitLevel(z_temp(i+1:i+round((1/(folder-1))*i)),Z(i+1:i+round((1/(folder-1))*i)));%测试模型的拟合度    R_train(step)=r_train;
R_test(step)=r_test;
r_total=weight*r_train+(1-weight)*r_test;    %整体拟合度
R_total(step)=r_total;
Total(step)=total;        %绘制R_total随样本数变化的动态图
set(h1,'xData',[1:step],'yData',R_total);
set(h2,'xData',[1:step],'yData',Total);
drawnow;
pause(0.1)
if r_max<r_total          %寻最优样本点
r_max=r_total;
Result(1,1)=i;
Result(2,1)=r_train;
Result(3,1)=r_test;
end
if r_train>train & r_test>test    %寻理想样本点列
Result(1,r_col)=i;
Result(2,r_col)=r_train;
Result(3,r_col)=r_test;
r_col=r_col+1;
end
step=step+1
end
%绘制返回结果图形
subplot(1,3,1)  %训练拟合度—样本数
plot(I,R_train)
hold on
plot([I(1),I(length(I))],[train,train])
axis([I(1),I(length(I)),-inf,inf])
title 'R train - number'
subplot(1,3,2)    %测试拟合度—样本数
plot(I,R_test)
hold on
plot([I(1),I(length(I))],[test,test])
axis([I(1),I(length(I)),-inf,inf])
title 'R test - number'
subplot(1,3,3)    %整体拟合度—样本数
plot(I,R_total)
hold on
plot([I(1),I(length(I))],[total,total])
axis([I(1),I(length(I)),-inf,inf])
title 'R total - number'
format bank
%以下是⼏点程序说明
%（1）次函数运⾏还伴随着整体拟合度r_total随样本数增⼤⽽变化的动态效果，把它注释掉也不会影响程序运⾏结果。%（2）我选拟合度指标不⼀定是最优的，如果还有更好的指标，可在FitLevel函数中修改。
%（3）交叉验证是⼀种验证模型精度的思想，不仅可⽤于检测插值模型精度，还可⽤于回归、分类、判别等模型中。%（4）⽤K-1份数据建模，对第K份进⾏验证会导致最后K份数据⽆法加⼊最优样本数的选举中，可以通过增加折数
%    和减⼩步长改进，不过这样会减缓程序速度。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。