一个简单的例子学明白用Python插值

这篇文章尝试通过一个简单的例子来为读者讲明白怎样使用Python实现数据插值。总共分3部分来介绍：

1. 为什么需要做插值这种事？
2. 通过拉格朗日插值法来看看插值这个事的理论要怎么理解？
3. Python实现拉格朗日插值的一个例子。

为什么需要做插值这种事？

这个答案很简单，无非两条：

首先，这个点上它没有数据或者数据不能用：（1）没采集到这个点上的数据；（2）采集到这个点上的数据了，但是数据明显是错误的。

其次，如果这个点上它没有数据的话呢，会对我们的建立的数据模型产生不好的影响，我们不得不想办法在这个缺失的点上给它想办法插上一个数据。

这个事，有点像下面这个桶桶，它少一个板就没法用，那我们想办法再找一个板板给它插到原来的地方，虽然可能没有原装的板板那么美观好看，但好歹这个桶桶它能用了呀。

那我们又不是神仙，怎样知道不存在数据的那个点上它的数据是个什么鬼呢？好在世间万物都不是孤立存在的，都会存在或多或少的联系，说的直白一点我们就靠这些联系来把这个缺少的值给猜出来。

还是以上面那个桶为例，我们要插入的多高多长的板子才能保证桶子盛的水最多且不漏出来呢？我们可以量一量缺少的那个板子它前后的板子的高度，然后取这两个高度的中间值作为插入的板子的高度，或者用整个桶的板子的平均高度，还有什么众数啊、固定值啊什么的能整的都可以整。

或者我们定义一个看上去比较NB的算法公式来确定这个板子的高度，比如用回归方法、拉格朗日插值法。那接下来我们一起看看拉格朗日插值，它其实也是一个非常简单的事。

通过拉格朗日插值法来看看插值这个事的理论要怎么理解？

我觉得用外国数学家的人名来命名一些方法、数学定理什么的是一件非常cd的事情，这极大的提高了我们中国学生的学习门槛，让很多简单的不能再简单的数学问题看上去那么高深莫测。好比我们说欧几里得空间，听着挺高大上的，让人望而生畏，其实它也就是种向量空间而已。

同样的，拉格朗日插值法实质上就是一种多项式插值法。而多项式插值法说的是，如果有n个点，每个点都有个x值对应的y值。那就必然能找到一个n-1次多项式使得这n个点的（x，y）值代入这个多项式成立。最简单的，好比说平面坐标上的两个点，必然能找的一个1次的式子y=kx+b满足这两个点的坐标值，更直白一点说，平面坐标上两点决定了一条直线。

那么，如果我们手上有n+1个点，但是只有其中n个点的（x，y）值，有一个点呢只知道它的x值，不知道它的y值。如果我们能找到满足已知的n个点的n-1次多项式（即插值公式），那我们就可以用这个多项式来算剩下的一个点的y值了。

拉格朗日先生呢，找到一种操作性比较强的办法确定这个多项式，所以我们把这个办法叫做拉格朗日插值法。

关于这个公式确定办法我会另外写文章讲，这里不再重复。

而且，对于实际工程应用来说，其实只要能够知道插值的思想，然后知道插值函数的调用方法即可，不需要记这些复杂的公式，毕竟我们主要目标是解决问题而不是去做数学家。

下面通过一个例子来说明Python进行数据插值的一般步骤。

Python实现拉格朗日插值的一个例子。

我们以后面参考资料中的一组数据为例来说明，需要数据源的朋友可以留言或私信我。

这组数据呢，是一个餐厅某段时间内的销量情况。数据源在excel中，我们使用pandas的read_excel方法将它读出来，放到一个dataframe中。实现代码如下：

import pandas as pd
from scipy.interpolate import lagrange #导入拉格朗日插值函数
import matplotlib.pyplot as plt

"""
读Excel文件
"""
inputfile = '../data/catering_sale.xls'#注意数据源excel文件的存放路径
data = pd.read_excel(inputfile)

然后呢，我们粗暴地假设这些销量中小于400或者大于5000的都是异常值，把它们置为空。

data[u'销量'][(data[u'销量']<400)|(data[u'销量']>5000)]=None#过滤异常值，将其变为空None
data_original = data.copy()#我们把原始数据保存一下

当然这个上下限400、5000，我们可以使用箱型图法等一些数据分布分析的办法来确定。

下面的代码就是对缺省的值进行插值了，你看其实代码很少。我们来分析一下。

前面这个ployinterp_column函数，是我们定义的插值函数，真正插值的操作是在后面的那两个嵌套的for循环中。

#定义列向量的插值函数
def ployinterp_column(s, n, k=4):
    y = s[list(range(n-k,n))+list(range(n+1,n+1+k))]#取一部分输入参数的值
    y = y[y.notnull()] #空值剔除掉
    return lagrange(y.index, list(y))(n)

for i in data.columns:
    for j in range(len(data)):
        if(data[i].isnull())[j]:  #如果为空
            data[i][j] = ployinterp_column(data[i],j)

上面代码外层的for循环中，for i in data.columns，是对data的每一列进行遍历，实质上我们这里只有两列，而且我们只需要对销量数据进行插值，所以我们这层可以省略。改成下面这个样子效果也是一样的：

for j in range(len(data)):
    if(data[u'销量'].isnull())[j]:  #如果为空
        data[u'销量'][j] = ployinterp_column(data[u'销量'],j)

我们遍历“销量”这一列，如果某个单元上为空，我们就调用ployinterp_column函数来对它进行插值。

注意到这个插值函数有3个参数，一个是我们要插值的整个列s，另一个是这列中为空的那个单元格的坐标n，还有一个k是我们取的整列中控制坐标n附近的几个值来进行插值（这里默认为4）。

插值前后的dataframe的比较如下图所示，我们在原来nan的位置上都自动的插入了一个值，而且这个值看上去还挺像那么回事的。

插值前后的对比

python里面实现拉格朗日插值很简单，直接调用scipy.interpolate里面的lagrange函数即可，但是需要注意的是我们在ployinterp_column函数中对k的取值的选择，k不宜取值过大或过小，过大的话会出现所谓的龙格现象。如下面两个图所示，k分别取4和5之后插值的效果，取5时有一个值时-70000多，明显是一个错误的。

k取4时的插值结果

k取5时的插值结果

所以，k的取值会影响插值的效果，而k具体取什么值合适，一般都是通过经验反复尝试几次来确定。

参考资料：

张良均等著，《Python与数据挖掘实践》