Python SciPy庫——擬合與插值

1.最小二乘擬合

實例1

import numpy as npnimport matplotlib.pyplot as pltnfrom scipy.optimize import leastsqnnnplt.figure(figsize=(9,9))nx=np.linspace(0,10,1000)nX = np.array([8.19, 2.72, 6.39, 8.71, 4.7, 2.66, 3.78])nY = np.array([7.01, 2.78, 6.47, 6.71, 4.1, 4.23, 4.05])n#計算以p為參數的直線和原始數據之間的誤差ndef f(p):n k, b = pn return(Y-(k*X+b))n#leastsq使得f的輸出數組的平方和最小，參數初始值為[1,0]nr = leastsq(f, [1,0])nk, b = r[0]nprint("k=",k,"b=",b)nnplt.scatter(X,Y, s=100, alpha=1.0, marker=o,label=u數據點)nny=k*x+bnnax = plt.gca()nnax.set_xlabel(..., fontsize=20)nax.set_ylabel(..., fontsize=20)n#設置坐標軸標籤字體大小nnplt.plot(x, y, color=r,linewidth=5, linestyle=":",markersize=20, label=u擬合曲線)nnplt.legend(loc=0, numpoints=1)nleg = plt.gca().get_legend()nltext = leg.get_texts()nplt.setp(ltext, fontsize=xx-large)nnplt.xlabel(u安培/A)nplt.ylabel(u伏特/V)nnplt.xlim(0, x.max() * 1.1)nplt.ylim(0, y.max() * 1.1)nnplt.xticks(fontsize=20)nplt.yticks(fontsize=20)n#刻度字體大小nnnplt.legend(loc=upper left)nnplt.show()n

實例2

#最小二乘擬合實例nimport numpy as npnfrom scipy.optimize import leastsqnimport pylab as plnndef func(x, p):n """n 數據擬合所用的函數: A*cos(2*pi*k*x + theta)n """n A, k, theta = pn return A*np.sin(k*x+theta) nndef residuals(p, y, x):n """n 實驗數據x, y和擬合函數之間的差，p為擬合需要找到的係數n """n return y - func(x, p)nnx = np.linspace(0, 20, 100)nA, k, theta = 10, 3, 6 # 真實數據的函數參數ny0 = func(x, [A, k, theta]) # 真實數據ny1 = y0 + 2 * np.random.randn(len(x)) # 加入雜訊之後的實驗數據 nnp0 = [10, 0.2, 0] # 第一次猜測的函數擬合參數nn# 調用leastsq進行數據擬合n# residuals為計算誤差的函數n# p0為擬合參數的初始值n# args為需要擬合的實驗數據nplsq = leastsq(residuals, p0, args=(y1, x))nnprint (u"真實參數:", [A, k, theta] )nprint (u"擬合參數", plsq[0]) # 實驗數據擬合後的參數nnpl.plot(x, y0, color=r,label=u"真實數據")npl.plot(x, y1, color=b,label=u"帶雜訊的實驗數據")npl.plot(x, func(x, plsq[0]), color=g, label=u"擬合數據")npl.legend()npl.show()n

2. 插值

# -*- coding: utf-8 -*-n"""nCreated on Thu Jul 27 16:42:30 2017nn@author: Delln"""nimport numpy as npnimport pylab as plnfrom scipy import interpolate nimport matplotlib.pyplot as pltnnx = np.linspace(0, 2*np.pi+np.pi/4, 10)ny = np.sin(x)nnx_new = np.linspace(0, 2*np.pi+np.pi/4, 100)nf_linear = interpolate.interp1d(x, y)ntck = interpolate.splrep(x, y)ny_bspline = interpolate.splev(x_new, tck)nnplt.xlabel(u安培/A)nplt.ylabel(u伏特/V)nnplt.plot(x, y, "o", label=u"原始數據")nplt.plot(x_new, f_linear(x_new), label=u"線性插值")nplt.plot(x_new, y_bspline, label=u"B-spline插值")nnpl.legend()npl.show()n

實例分析

# -*- coding: utf-8 -*-n"""nCreated on Thu Jul 27 16:53:21 2017nn@author: Delln"""nnimport numpy as npnfrom scipy import interpolatenimport pylab as pln#創建數據點集並繪製npl.figure(figsize=(12,9))nx = np.linspace(0, 10, 11)ny = np.sin(x)nax=pl.plot()nnpl.plot(x,y,ro)n#建立插值數據點nxnew = np.linspace(0, 10, 101)nfor kind in [nearest, zero,linear,quadratic]:n #根據kind創建插值對象interp1dn f = interpolate.interp1d(x, y, kind = kind)n ynew = f(xnew)#計算插值結果n pl.plot(xnew, ynew, label = str(kind))nnpl.xticks(fontsize=20)npl.yticks(fontsize=20)nnpl.legend(loc = lower right)npl.show()n

B樣條曲線插值

一維數據的插值運算可以通過 interp1d()實現。

其調用形式為：

Interp1d可以計算x的取值範圍之內任意點的函數值，並返回新的數組。

interp1d(x, y, kind=『linear』, …)

參數 x和y是一系列已知的數據點

參數kind是插值類型，可以是字元串或整數

B樣條曲線插值

Kind給出了B樣條曲線的階數：

? 『

zero『 『nearest』 :0階梯插值，相當於0階B樣條曲線

? 『slinear』『linear』 :線性插值，相當於1階B樣條曲線

? 『quadratic』『cubic』:2階和3階B樣條曲線，更高階的曲線可以直接使用整數值來指定

（1）#創建數據點集：

import numpy as np

x = np.linspace(0, 10, 11)

y = np.sin(x)

（2）#繪製數據點集：

import pylab as pl

pl.plot(x,y,ro)

創建interp1d對象f、計算插值結果：

xnew = np.linspace(0, 10, 11)

from scipy import interpolate

f = interpolate.interp1d(x, y, kind = kind)

ynew = f(xnew)

根據kind類型創建interp1d對象f、計算並繪製插值結果：

xnew = np.linspace(0, 10, 11)

for kind in [nearest, zero,linear,quadratic]:

#根據kind創建插值對象interp1d

f = interpolate.interp1d(x, y, kind = kind)

ynew = f(xnew)#計算插值結果

pl.plot(xnew, ynew, label = str(kind))#繪製插值結果

如果我們將代碼稍作修改增加一個5階插值

# -*- coding: utf-8 -*-n"""nCreated on Thu Jul 27 16:53:21 2017nn@author: Delln"""nnimport numpy as npnfrom scipy import interpolatenimport pylab as pln#創建數據點集並繪製npl.figure(figsize=(12,9))nx = np.linspace(0, 10, 11)ny = np.sin(x)nax=pl.plot()nnpl.plot(x,y,ro)n#建立插值數據點nxnew = np.linspace(0, 10, 101)nfor kind in [nearest, zero,linear,quadratic,5]:n #根據kind創建插值對象interp1dn f = interpolate.interp1d(x, y, kind = kind)n ynew = f(xnew)#計算插值結果n pl.plot(xnew, ynew, label = str(kind))nnpl.xticks(fontsize=20)npl.yticks(fontsize=20)nnpl.legend(loc = lower right)npl.show()n運行得到n

發現5階已經很接近正弦曲線，但是如果x值選取範圍較大，則會出現跳躍。

關於擬合與插值的數學基礎可參見霍開拓：擬合與插值的區別？

左邊插值，右邊擬合

仔細看有啥不一樣

插值曲線要過數據點，擬合曲線整體效果更好。

插值，對準了才可以插嗎，那就一定得過數據點。擬合，就是要得到最接近的結果，是要看總體效果。

既然理想(思路)不一樣，那麼三觀和行為（特點和策略）也就不一樣啦。

插值是指已知某函數的在若干離散點上的函數值或者導數信息，通過求解該函數中待定形式的插值函數以及待定係數，使得該函數在給定離散點上滿足約束。

所謂擬合是指已知某函數的若干離散函數值{f1,f2,…,fn}，通過調整該函數中若干待定係數f(λ1, λ2,…,λn), 使得該函數與已知點集的差別(最小二乘意義)最小。如果待定函數是線性，就叫線性擬合或者線性回歸(主要在統計中)，否則叫作非線性擬合或者非線性回歸。表達式也可以是分段函數，這種情況下叫作樣條擬合。

從幾何意義上將，擬合是給定了空間中的一些點，找到一個已知形式未知參數的連續曲面來最大限度地逼近這些點；而插值是找到一個( 或幾個分片光滑的)連續曲面來穿過這些點。

插值法

以下引自某科

Lagrange插值

Lagrange插值是n次多項式插值，其成功地用構造插值基函數的 方法解決了求n次多項式插值函數問題。

基本思想　將待求的n次多項式插值函數pn(x）改寫成另一種表示方式，再利用插值條件⑴確定其中的待定函數，從而求出插值多項式。

Newton插值

Newton插值也是n次多項式插值，它提出另一種構造插值多項式的方法，與Lagrange插值相比，具有承襲性和易於變動節點的特點。

基本思想　將待求的n次插值多項式Pn（x）改寫為具有承襲性的形式，然後利用插值條件⑴確定Pn（x）的待定係數，以求出所要的插值函數。

Hermite插值

Hermite插值是利用未知函數f(x)在插值節點上的函數值及導數值來構造插值多項式的，其提法為：給定n+1個互異的節點x0,x1，……,xn上的函數值和導數值求一個2n+1次多項式H2n+1(x)滿足插值條件H2n+1(xk)=ykH2n+1(xk)=yk k=0,1,2，……，n ⒀如上求出的H2n+1(x）稱為2n+1次Hermite插值函數，它與被插函數一般有更好的密合度.

基本思想利用Lagrange插值函數的構造方法，先設定函數形式，再利用插值條件⒀求出插值函數.

貌似插值節點取的越多，差值曲線或曲面越接近原始曲線/曲面，因為採樣多嘛。但事實總是不像廣大人民群眾想的那樣，隨著插值節點的增多，多項式次數也在增高，插值曲線在一些區域出現跳躍，並且越來越偏離原始曲線。這個現象被 Tolmé Runge 發現並解釋，然後就以他的名字命名這種現象。It was discovered by Carl David Tolmé Runge (1901) when exploring the behavior of errors when using polynomial interpolation to approximate certain functions.

為了解決這個問題，人們發明了分段插值法。分段插值一般不會使用四次以上的多項式，而二次多項式會出現尖點，也是有問題的。所以就剩下線性和三次插值，最後使用最多的還是線性分段插值，這個好處是顯而易見的。

擬合

最小二乘

如何找到最接近原始曲線或者數據點的擬合曲線，這不是一件容易操作的事。要想整體最接近，直接的想法就是擬合曲線的每一點到原始曲線的對應點的最接近，簡單點說就是兩曲線上所有點的函數值之差的絕對值之和最小。看似解決問題，但絕對值在數學上向來是個不好交流的語言障礙患者，那然後又該怎麼辦。數學家說了既然辦不了你絕對值之和，那就辦了你家親戚，就看你平方之和長得像。於是就找了這個長得像的來背黑鍋，大家都表示很和諧。然後給這種操作冠之名曰最小二乘法。

官方一點的表述 ， 選擇參數c使得擬合模型與實際觀測值在曲線擬合各點的殘差(或離差)ek=yk－f(xk，c)的加權平方和達到最小，此時所求曲線稱作在加權最小二乘意義下對數據的擬合曲線，這種方法叫做最小二乘法。

推薦閱讀：