這篇文章嘗試通過一個簡單的例子來為讀者講明白怎樣使用Python實現數據插值。總共分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中。實現代碼如下:
importpandasaspdfromscipy.interpolateimportlagrange#導入拉格朗日插值函數importmatplotlib.pyplotasplt"""讀Excel文件"""inputfile='../data/catering_sale.xls'#注意數據源excel文件的存放路徑data=pd.read_excel(inputfile)然後呢,我們粗暴地假設這些銷量中小於400或者大於5000的都是異常值,把它們置為空。
data[u'銷量'][(data[u'銷量']<400)|(data[u'銷量']>5000)]=None#過濾異常值,將其變為空Nonedata_original=data.copy()#我們把原始數據保存一下當然這個上下限400、5000,我們可以使用箱型圖法等一些數據分布分析的辦法來確定。
下面的代碼就是對預設的值進行插值了,你看其實代碼很少。我們來分析一下。
前面這個ployinterp_column函數,是我們定義的插值函數,真正插值的操作是在後面的那兩個嵌套的for循環中。
#定義列向量的插值函數defployinterp_column(s,n,k=4):y=s[list(range(n-k,n))+list(range(n+1,n+1+k))]#取一部分輸入參數的值y=y[y.notnull()]#空值剔除掉returnlagrange(y.index,list(y))(n)foriindata.columns:forjinrange(len(data)):if(data[i].isnull())[j]:#如果為空data[i][j]=ployinterp_column(data[i],j)上面代碼外層的for循環中,for i in data.columns,是對data的每一列進行遍歷,實質上我們這裡只有兩列,而且我們只需要對銷量數據進行插值,所以我們這層可以省略。改成下面這個樣子效果也是一樣的:
forjinrange(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與數據挖掘實踐》