未来可以预测吗？（Matlab系列之插值与拟合)

前言

在前几天看了一个科幻的影视剧，名字叫“开发者”，有8集，时间也不短，还没来得及看完，但是在开始的几集里头，就出现了一种技术。

大概是这样的：影视里有个程序猿，他写出了一种程序，可以预测米线虫未来几十秒内的运动，这一部分片段出现的较多的词是同步以及相关性，不过按我猜测，这个预测的功能就是通过拟合来实现的，这也是为什么要用米线虫来做实验，因为是种简单的有机生物体，不会像人一样有太多的想法，此外，影视里还有将过去发生的事进行再次展现的片段等等，有兴趣的朋友可以去看看。

友情提示：开场有点惊悚，想看但又有点害怕的，就找上你的小伙伴一起看吧~

介绍

扯了一些东西，开始今天的内容吧，插值与拟合是数据分析中常见的方法，在解决实际问题的工程实践和科学实验过程中，经常会需要通过研究某些变量之间的函数关系来让我们了解事物的内在规律以及本质属性，这些变量之间的存在的未知的函数关系通常又隐藏于试验以及观测的数据中，因此通过数据找出函数关系是解决问题的关键，开始说的那个电影片段就可以这样举例，程序猿通过寻找米线虫的运动规律，从而预测了未来几十秒米线虫的运动轨迹。

MATLAB中常用的插值和拟合函数如下表所示：

函数	描述
interp1(x,y,xo)	1维线性插值
interp1(x,y,xo,‘spline’)	1维3次样条插值
interp1(x,y,xo,‘cubic’)	1维3次插值
interp2(x,y,Z,xi,yi)	2维线性插值
interp2(x,y,Z,xi,yi,‘cubic’)	2维3次插值
interp2(x,y,Z,xi,yi,‘nearest’)	2维最近邻插值
polyfit(x,y,n)	对给定的数据进行n次多项式拟合

插值

一维插值

先简单解释下一维插值在数学中的描述，假设已知一个函数的表达式是很复杂的或者是未知的，在一个区间内，n个变量对应n个函数值（就相当于xy的坐标轴，x在y上有一个对应的值，即呈现出函数关系），且该函数在这个区间内有若干阶次的导数，求出在该区间内，x轴上任意一点对应的y轴值（这是个近似值）。

表述大概就是上面那个样子，如果还不太理解，就需要各位自己找相关资料深入了解了，接下来直接看MATLAB的操作。

y=interp1(x0,y0,x,'method')

注：method是指定的插值方法，如果没填，默认使用线性插值“linear”，其他的方法有最近项插值“nearest”，立方样条插值“spline”，立方插值“cubic”，要注意的一点，x0需要是单调的，假如x0还是等距的数值，则可以使用快速插值法，快速插值就是在以上的插值法前加一个“*”，例如：

y=interp1(x0,y0,x,'*nearest')

接下来上程序演练一下：

x=0:0.1:10;
y=sin(x);%产生函数关系的数据
xi=0:0.25:10;%需要插值的点
yi=interp1(x,y,xi);%使用插值函数得出插值点对应的值
plot(x,y,'k*',xi,yi,'-r');%画出原图与插值后的图进行对比

结果，黑色的“*”就是原图，红色的连续线条就是插值结果：

二维插值

二维插值和一维的概念相似，只是将一维转变成了二维的形式，原本是一段区间，二维则是一块矩形区域，数学表述：已知二元函数g(x,y)在某矩形区域R[a,b]x[c,d]上互异节点(x[i],y[j])的函数值为z[i,j]，求出区域R内任一点对应的函数值g(x,y)的近似值。

可以用函数式这样理解：

一维：y=f(x)；注：x∈[a,b]
二维：z(x,y)=g(x,y); 注：(x,y ∈ R[a,b]x[c,d])

函数使用：

z=interp2(x0,y0,z0,x,y,‘method’)

其中，x0、y0分别是m维和n维的向量，表示节点；z0则是一个n x m维的矩阵，表示节点值；x、y是一维数组，表示插值点，两者为不同方向的向量，比如：x个是行向量，y则是列向量；z则是得到的插值，是一个矩阵；method和一维的一致。

来个程序看看：

[X,Y,V] = peaks(10); %产生10*10的矩阵
[Xq,Yq] = meshgrid(-3:.1:3,-3:.1:3);%产生新的向量组：X和Y
Vq = interp2(X,Y,V,Xq,Yq); %对新的向量组进行二维插值
%绘图做对比
figure;
subplot(211);mesh(X,Y,V);%绘制原始矩阵的立体图
subplot(212);mesh(Xq,Yq,Vq);%绘制二维插值结果的立体图

结果：

可以看到插值成功的将区域内非X和Y的部分也成功计算出，数据量明显增多，如果看图看不明白，那看下变量的数据量也可得知：

拟合

具有“预测”功能的东西来了，还是先简单进行下数学描述：

已知一组二维的数据（xi,yi）,其中i=1，2，3，…，n，xi各不相同，求得一个函数f(x)，使得在某种情况下该函数会和所有的数据点接近，即曲线拟合。

在MATLAB中，函数polyfit可以解决最小二乘曲线拟合问题，能对给定数据进行多项式拟合，最好给出多项式的系数，使用方法如：

p=polyfit(x,y,n)

注：采用n次多项式p对数据x和y进行拟合，从而使得p(x)与y的最小均方差最小，x和y为给定的数据，n为多项式拟合的次数，p就是拟合出的多项式系数。

多项式就可以根据系数进行表示：

y=pm*xm+p(m-1)*x(m-1)...+p1*x+p0;%系数p就等于{pm,pm-1,....,p1,p0}

想要求得在x在拟合多项式上的对应值y，则可用polyval进行计算。

程序举例：

%步骤1 产生拟合所需数据并拟合出对应多项式
x=(0:0.1:2.5)';
y=erf(x);%就一个error function，
p=polyfit(x,y,6)%拟合6次，求得对应的多项式系数
%步骤2 用同样函数产生的数据结果与拟合多项式计算的结果进行对比
xi=(0:0.1:5)';
yi=erf(xi);%使用同一个函数构建一个新的数据用以做对比
s=polyval(p,xi);%使用拟合多项式计算函数值
plot(xi,yi,'b*',xi,s,'-r');%绘制对比结果图
axis([0 5 0 2]);%限定绘图区间

结果：

蓝色的就是原函数得到的结果，红色的就是通过0_{2.5区间的数据得到的拟合多项式，再将该多项式应用0}5的区间的结果，可以明显看到拟合的结果在3之前基本一致，在4之后才完全偏离原函数的结果，很显然“预测”成功。

所以，如果数据量大一些~

是不是真的可以预测未来了呢？

实例应用

假设某个小企业的所得利润如下：

年份	2001	2002	2003	2004	2005	2006	2007
利润（万元）	70	123	145	155	176	196	211

问：预测2008年和2009年的利润

程序如下：

x=[2001 2002 2003 2004 2005 2006 2007];
y=[70 123 145 155 176 196 211];
p=polyfit(x,y,1);
xn=[2008 2009];
yn=polyval(p,xn)

不知道有没有发现，这是个单调上升的盈利企业…

好了，本篇记录到此结束~

---

---

更多精彩，等你发现~

---

---