在離散數據的基礎上補插連續函數，使得這條連續曲線通過全部給定的離散數據點。插值是離散函數逼近的重要方法，利用它可通過函數在有限個點處的取值狀況，估算出函數在其他點處的近似值。曲面插值是對三維數據進行離散逼近的方法，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),首先用訓練集對分類器進行訓練,再利用驗證集來測試訓練得到的模型(model),以此來做為評價分類器的性能指標。在此基礎上，反復的做訓練、測試以及模型的選擇。包含簡單交叉驗證，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。
%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.com
%設置參數默認值
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份數據無法加入最優樣本數的選舉中，可以通過增加折數
% 和減小步長改進，不過這樣會減緩程序速度。