[Matlab数据分析]使用有限数据推算其他数值

2支棒棒糖

插值与拟合是使用有限数据对其他数值进行推算的基本方法。

插值可以通过有限点来建立简单连续的解析模型，并根据该模型得到未知点处的值。插值的方法很多，插值按对象可分为一维插值、多维插值和离散数据插值。

1．一维插值

MATLAB提供了interp1函数实现一维插值，该函数的调用格式如下：

1. yi = interp1(x,Y,xi)
   
2. yi = interp1(Y,xi)
   
3. yi = interp1(x,Y,xi,method)
   
4. yi = interp1(x,Y,xi,method,'extrap')
   
5. yi = interp1(x,Y,xi,method,extrapval)
   
6. pp = interp1(x,Y,method,'pp')

复制代码

其中，x为自变量的取值向量；Y为对应的函数值；xi为插值点；yi为插值结果；'extrap'表示进行外插，extrapval表示外插标量；'pp'为采用分段多项式进行插值，pp为得到的结果。

另外，MATLAB还提供了其他插值函数，如pchip、spline等；这里暂不对这些函数做具体介绍。

例10-6，一维插值示例。

在命令行窗口输入：

1. x = 0:10;
   
2. y = sin(x);
   
3. xi = 0:.25:10;
   
4. y1 = interp1(x,y,xi);                              %线性插值
   
5. subplot(121);plot(x,y,'o',xi,y1);title('linear')
   
6. y2 = interp1(x,y,xi,'spline');                     %三次样条插值
   
7. subplot(122);plot(x,y,'o',xi,y2) ;title('spline')

复制代码

MATLAB还提供了griddedInterpolant函数用于栅格的一维插值，该函数的调用格式如下所示：

1. F = griddedInterpolant(x,v)
   
2. F = griddedInterpolant(X1,X2,...,Xn,V)
   
3. F = griddedInterpolant(V)
   
4. F = griddedInterpolant({xg1,xg2,...,xgn},V)
   
5. F = griddedInterpolant(...,Method)

复制代码

其中，F为返回的插值函数；method为插值方法，如表10.1所示；x、v为同大小向量输入；X1, X2, ..., Xn构成n维数组，构成插值栅格点；V与X1, X2, ..., Xn构成的数组同维，为栅格点的值；{xg1,xg2,...,xgn}为单元数组。

例10-7，多列一维插值示例。

在命令行窗口输入：

1. x = (1:5)';
   
2. V = [x, 2*x, 3*x]
   
3. samplePoints = {x, 1:size(V,2)};
   
4. F = griddedInterpolant(samplePoints,V);
   
5. queryPoints = {(1:0.5:5),1:size(V,2)};
   
6. Vq = F(queryPoints)

复制代码

输出结果如下：

1. V = 1       2    3
   
2.        2     4     6
   
3.        3     6     9
   
4.        4     8    12
   
5.        5    10     15
   
6. Vq =1.0000     2.0000     3.0000
   
7.          1.5000     3.0000     4.5000
   
8.          2.0000     4.0000     6.0000
   
9.          2.5000     5.0000     7.5000
   
10.          3.0000     6.0000     9.0000
    
11.          3.5000     7.0000    10.5000
    
12.          4.0000     8.0000    12.0000
    
13.          4.5000     9.0000    13.5000
    
14.          5.0000    10.0000    15.0000

复制代码

提示:上例中实际进行了二维插值，然后再从二维插值中取样得到一维插值。对于多列向量，这种插值方法是一种高效的方法。

2．多维插值

MATLAB提供多维插值的函数包括griddedInterpolant、interp2、interp3、interpn等。下面介绍使用这些函数的插值方法。

前文中已经提到的griddedInterpolant可以进行二维插值，如例10-6所示。

例10-8，细化网格采样点示例。

在命令行窗口输入：

1. [X1, X2] = ndgrid(1:10,1:10);
   
2. V = X1.^2 + X2.^2;
   
3. subplot(121);mesh(X1,X2,V);
   
4. F = griddedInterpolant(X1,X2,V, 'cubic');
   
5. [Xq1, Xq2] = ndgrid(1:0.25:10,1:0.25:10);
   
6. Vq = F(Xq1,Xq2);
   
7. subplot(122);mesh(Xq1,Xq2,Vq);  

复制代码

MATLAB提供了interp2函数用于二维数据的插值，该函数的调用格式如下所示：

1. ZI = interp2(X,Y,Z,XI,YI)
   
2. ZI = interp2(Z,XI,YI)
   
3. ZI = interp2(Z,ntimes)
   
4. ZI = interp2(X,Y,Z,XI,YI,method)
   
5. ZI = interp2(...,method, extrapval)

复制代码

其中，X、Y、Z为采样数据点坐标；XI、YI、ZI为待插值点的坐标和返回的值；method为插值方法，见上文表10.1部分说明；extrapval表示外插标量。

例10-9，使用interp2进行二维插值示例。

在命令行窗口输入：

1. [X,Y] = meshgrid(-3:.5:3);
   
2. Z = peaks(X,Y);
   
3. [XI,YI] = meshgrid(-3:.25:3);
   
4. ZI1 = interp2(X,Y,Z,XI,YI);
   
5. ZI2 = interp2(X,Y,Z,XI,YI, 'nearest');
   
6. ZI3 = interp2(X,Y,Z,XI,YI,'spline');
   
7. subplot(221);surf(X,Y,Z), axis([-3 3 -3 3 -5 5]) ;title('Sample')      
   
8. subplot(222);surf(XI,YI,ZI1), axis([-3 3 -3 3 -5 5]) ;title('Linear')        
   
9. subplot(223);surf(XI,YI,ZI2), axis([-3 3 -3 3 -5 5]) ;title('Nearest')      
   
10. subplot(224);surf(XI,YI,ZI3), axis([-3 3 -3 3 -5 5]) ;title('Spline')

复制代码

MATLAB提供了interp3函数用于三维数据的插值，该函数的调用格式如下所示：

1. VI = interp3(X,Y,Z,V,XI,YI,ZI)
   
2. VI = interp3(V,XI,YI,ZI)
   
3. VI = interp3(V,ntimes)
   
4. VI = interp3(...,method)
   
5. VI = interp3(...,method,extrapval)

复制代码

其中，X、Y、Z、V为采样数据点坐标和值；XI、YI、ZI、VI为待插值点的坐标和返回的值；method为插值方法，见上文中表10.1部分说明；extrapval表示外插标量；X、Y、Z默认为数据下标；ntimes为对数据点间插值的次数。

例10-10，三维插值示例。

在命令行窗口输入：

1. [x,y,z,v] = flow(10);
   
2. subplot(121);slice(x,y,z,v,[6 9.5],2,[-2 .2])
   
3. [xi,yi,zi] = meshgrid(.1:.25:10, -3:.25:3, -3:.25:3);
   
4. vi = interp3(x,y,z,v,xi,yi,zi,'spline'); % vi is 25-by-40-by-25
   
5. subplot(122);slice(xi,yi,zi,vi,[6 9.5],2,[-2 .2])

复制代码

MATLAB提供了interpn函数用于多维数据的插值，由于前文提供的函数已基本能满足使用要求，故本书不再介绍此函数。如有需要，可参考帮助文档。

3．离散数据插值

离散数据是点X与值V对应的非连续数据。由于离散数据的点没有规律可言，所以相对而言插值较为复杂。

MATLAB提供的离散数据插值函数包括griddata、griddatan和TriScatteredInterp函数。griddata函数的调用格式如下：

1. vq = griddata(x,y,v,xq,yq)
   
2. vq = griddata(x,y,z,v,xq,yq,zq)
   
3. vq = griddata(..., method)

复制代码

其中，x、y、z为坐标，v为对应采样值；xq、yq、zq为插值点坐标，vq为插值结果；method为插值方法，可选'linear'、'cubic'、'natural'、'nearst'和'v4'，'natural'为自然相邻插值，'v4'为MATLAB4采用的方法，其余可参考前文同样的参数。

例10-11，使用griddata函数进行离散数据插值示例。

在命令行窗口输入：

1. xy = -2.5 + 5*gallery('uniformdata',[200 2],0);
   
2. x = xy(:,1); y = xy(:,2);
   
3. v = x.*exp(-x.^2-y.^2);
   
4. [xq,yq] = meshgrid(-2:.2:2, -2:.2:2);
   
5. vq = griddata(x,y,v,xq,yq); h = gca; set(h,'XLim',[-2.7 2.7]); set(h,'YLim',[-2.7   2.7]);
   
6. subplot(121); plot3(x,y,v,'o');
   
7. h = gca;set(h,'XLim',[-2.7 2.7]);set(h,'YLim',[-2.7 2.7]);
   
8. subplot(122);mesh(xq,yq,vq);
   
9. hold on;plot3(x,y,v,'o');
   
10. h = gca;set(h,'XLim',[-2.7 2.7]);set(h,'YLim',[-2.7 2.7]);

复制代码

TriScatteredInterp函数的调用格式为：

1. F = TriScatteredInterp()
   
2. F = TriScatteredInterp(X, V)
   
3. F = TriScatteredInterp(X, Y, V)
   
4. F= TriScatteredInterp(X, Y, Z, V)
   
5. F = TriScatteredInterp(..., method)

复制代码

其中，X、Y、Z为坐标，V为坐标对应值，F为返回值；method为插值方法，可选'linear'、'natural'或'nearst'。

例10-12，使用TriScatteredInterp函数进行离散数据插值示例。

在命令行窗口输入：

1. x = rand(100,1)*4-2;
   
2. y = rand(100,1)*4-2;
   
3. z = x.*exp(-x.^2-y.^2);
   
4. subplot(121);plot3(x,y,z,'o');  
   
5. F = TriScatteredInterp(x,y,z);
   
6. ti = -2:.25:2;
   
7. [qx,qy] = meshgrid(ti,ti);
   
8. qz = F(qx,qy);
   
9. subplot(122);mesh(qx,qy,qz);hold on;plot3(x,y,z,'o');

复制代码

4.曲线拟合

曲线拟合可以支持在数据间支持定量关系式，这在工程和科学应用上非常重要。本小节着重介绍在MATLAB中使用多项式进行拟合的方法。

MATLAB提供用于多项式拟合的函数为polyfit函数，该函数的调用格式如下：

p = polyfit(x,y,n)
其中，x、y、n为输入的x值、y值和拟合多项式的阶次；p为得到的多项式系数。

例10-13，使用polyfit函数进行曲线拟合示例。

在命令行窗口输入：

1. x = (0: 0.7: 5)';
   
2. x1= (0: 0.05: 5)';
   
3. y = sin(x);
   
4. p1 = polyfit(x,y,4)
   
5. p2 = polyfit(x,y,7)
   
6. p3 = polyfit(x,y,10)
   
7. f1 = polyval(p1,x1);
   
8. f2 = polyval(p2,x1);
   
9. f3 = polyval(p3,x1);
   
10. subplot(131);plot(x,y,'o',x1,f1,'-'); axis([0  5 -1.2 1.2])
    
11. subplot(132);plot(x,y,'o',x1,f2,'-'); axis([0  5 -1.2 1.2])
    
12. subplot(133);plot(x,y,'o',x1,f3,'-'); axis([0  5 -1.2 1.2])

复制代码

输出结果如下：

1. p1 =     0.0196     -0.1059     -0.2303      1.1746     -0.0080
   
2. p2 =     0.0001     -0.0030      0.0211     -0.0303     -0.1264    -0.0275  1.0073    -0.0000
   
3. p3 = Columns 1 through 9
   
4.         0.0002     -0.0030      0.0217     -0.0749      0.1150      0  -0.1857    0     0
   
5. Columns 10 through 11
   
6.         0.9719      0.0000

复制代码