Python SciPy庫——插值與擬合

插值與擬合

原文鏈接：https://zhuanlan.zhihu.com/p/28149195

 

1.最小二乘擬合

實例1

 

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
## 設置字符集，防止中文亂碼
import matplotlib
matplotlib.rcParams['font.sans-serif']=[u'simHei']
matplotlib.rcParams['axes.unicode_minus']=False

plt.figure(figsize=(9,9))
x=np.linspace(0,10,1000)
X = np.array([8.19, 2.72, 6.39, 8.71, 4.7, 2.66, 3.78])
Y = np.array([7.01, 2.78, 6.47, 6.71, 4.1, 4.23, 4.05])
#計算以p為參數的直線和原始數據之間的誤差
def f(p):
    k, b = p
    return(Y-(k*X+b))
#leastsq使得f的輸出數組的平方和最小，參數初始值為[1,0]
r = leastsq(f, [1,0])
k, b = r[0]
print("k=",k,"b=",b)

plt.scatter(X,Y, s=100, alpha=1.0, marker='o',label=u'數據點')

y=k*x+b

ax = plt.gca()

ax.set_xlabel(..., fontsize=20)
ax.set_ylabel(..., fontsize=20)
#設置坐標軸標簽字體大小

plt.plot(x, y, color='r',linewidth=5, linestyle=":",markersize=20, label=u'擬合曲線')

plt.legend(loc=0, numpoints=1)
leg = plt.gca().get_legend()
ltext  = leg.get_texts()
plt.setp(ltext, fontsize='xx-large')

plt.xlabel(u'安培/A')
plt.ylabel(u'伏特/V')

plt.xlim(0, x.max() * 1.1)
plt.ylim(0, y.max() * 1.1)

plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
#刻度字體大小


plt.legend(loc='upper left')

plt.show()

 

k= 0.6134953491930442 b= 1.794092543259387

 

實例2

 

# -*- coding: utf-8 -*-

#最小二乘擬合實例
import numpy as np
from scipy.optimize import leastsq
import pylab as pl
## 設置字符集，防止中文亂碼
import matplotlib
matplotlib.rcParams['font.sans-serif']=[u'simHei']
matplotlib.rcParams['axes.unicode_minus']=False

def func(x, p):
    """
    數據擬合所用的函數: A*cos(2*pi*k*x + theta)
    """
    A, k, theta = p
    return A*np.sin(k*x+theta)   

def residuals(p, y, x):
    """
    實驗數據x, y和擬合函數之間的差，p為擬合需要找到的系數
    """
    return y - func(x, p)

x = np.linspace(0, 20, 100)
A, k, theta = 10, 3, 6 # 真實數據的函數參數
y0 = func(x, [A, k, theta]) # 真實數據
y1 = y0 + 2 * np.random.randn(len(x)) # 加入噪聲之后的實驗數據    

p0 = [10, 0.2, 0] # 第一次猜測的函數擬合參數

# 調用leastsq進行數據擬合
# residuals為計算誤差的函數
# p0為擬合參數的初始值
# args為需要擬合的實驗數據
plsq = leastsq(residuals, p0, args=(y1, x))

print (u"真實參數:", [A, k, theta] )
print (u"擬合參數", plsq[0]) # 實驗數據擬合后的參數

pl.plot(x, y0, color='r',label=u"真實數據")
pl.plot(x, y1, color='b',label=u"帶噪聲的實驗數據")
pl.plot(x, func(x, plsq[0]), color='g', label=u"擬合數據")
pl.legend()
pl.show()

 

真實參數: [10, 3, 6]
擬合參數 [-1.16428658  0.24215786 -0.794681  ]

 

2.插值

實例1

 

# -*- coding: utf-8 -*-

# -*- coding: utf-8 -*-

import numpy as np
import pylab as pl
from scipy import interpolate 
import matplotlib.pyplot as plt
## 設置字符集，防止中文亂碼
import matplotlib
matplotlib.rcParams['font.sans-serif']=[u'simHei']
matplotlib.rcParams['axes.unicode_minus']=False

x = np.linspace(0, 2*np.pi+np.pi/4, 10)
y = np.sin(x)

x_new = np.linspace(0, 2*np.pi+np.pi/4, 100)
f_linear = interpolate.interp1d(x, y)
tck = interpolate.splrep(x, y)
y_bspline = interpolate.splev(x_new, tck)

plt.xlabel(u'安培/A')
plt.ylabel(u'伏特/V')

plt.plot(x, y, "o",  label=u"原始數據")
plt.plot(x_new, f_linear(x_new), label=u"線性插值")
plt.plot(x_new, y_bspline, label=u"B-spline插值")

pl.legend()
pl.show()

 

 

 

實例2

# -*- coding: utf-8 -*-

import numpy as np
from scipy import interpolate
import pylab as pl
## 設置字符集，防止中文亂碼
import matplotlib
matplotlib.rcParams['font.sans-serif']=[u'simHei']
matplotlib.rcParams['axes.unicode_minus']=False

#創建數據點集並繪制
pl.figure(figsize=(12,9))
x = np.linspace(0, 10, 11)
y = np.sin(x)
ax=pl.plot()

pl.plot(x,y,'ro')
#建立插值數據點
xnew = np.linspace(0, 10, 101)
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))

pl.xticks(fontsize=20)
pl.yticks(fontsize=20)

pl.legend(loc = 'lower right')
pl.show()

 

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階插值

import numpy as np
from scipy import interpolate
import pylab as pl
#創建數據點集並繪制
pl.figure(figsize=(12,9))
x = np.linspace(0, 10, 11)
y = np.sin(x)
ax=pl.plot()

pl.plot(x,y,'ro')
#建立插值數據點
xnew = np.linspace(0, 10, 101)
for kind in ['nearest', 'zero','linear','quadratic',5]:
    #根據kind創建插值對象interp1d
    f = interpolate.interp1d(x, y, kind = kind)
    ynew = f(xnew)#計算插值結果
    pl.plot(xnew, ynew, label = str(kind))

pl.xticks(fontsize=20)
pl.yticks(fontsize=20)

pl.legend(loc = 'lower right')
pl.show()
運行得到
 

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

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

左邊插值，右邊擬合