Python數(shù)學(xué)建模系列（四）：數(shù)值逼近
1. 一維插值
插值：求過(guò)已知有限個(gè)數(shù)據(jù)點(diǎn)的近似函數(shù)。
插值函數(shù)經(jīng)過(guò)樣本點(diǎn)，擬合函數(shù)一般基于最小二乘法盡量靠近所有樣本并穿過(guò)。常見(jiàn)差值方法有拉格朗日插值法、分段插值法、樣條插值法。
?interp1d(x, y) 計(jì)算一維插值
1.1 線性插值與樣條插值(B-spline)
例1：某電學(xué)元件的電壓數(shù)據(jù)記錄在0~2.25πA范圍與電流關(guān)系滿足正弦函數(shù)，分別用線性插值和樣條插值方法給出經(jīng)過(guò)數(shù)據(jù)點(diǎn)的數(shù)值逼近函數(shù)曲線。
Demo代碼
import matplotlib
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
# 引入中文字體
font = {
"family": "Microsoft YaHei"
}
matplotlib.rc("font", **font)
# 初始數(shù)據(jù)量 0 - 2.25pi 分為10份 均勻分
x = np.linspace(0, 2.25 * np.pi, 10)
y = np.sin(x)
# 得到差值函數(shù) （使用線性插值）
f_linear = interpolate.interp1d(x, y)
# 新數(shù)據(jù) 0 - 2.25pi 分為100份 均勻分 （線性插值）
x_new = np.linspace(0, 2.25 * np.pi, 100)
y_new = f_linear(x_new)
# 使用B-spline插值
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"原始數(shù)據(jù)")
plt.plot(x_new, f_linear(x_new), label=u"線性插值")
plt.plot(x_new, y_bspline, label=u"B-spline插值")
plt.legend()
plt.show()
輸出：
涉及知識(shí)點(diǎn)：
numpy.linspace scipy.interpolate.interp1d scipy.interpolate.splrep
1.2 高階樣條插值
隨著插值節(jié)點(diǎn)增多，多項(xiàng)式次數(shù)也變高，插值曲線在一些區(qū)域出現(xiàn)跳躍，并且越來(lái)越偏離原始曲線，稱為龍格現(xiàn)象。
例2：某電學(xué)元件的電壓數(shù)據(jù)記錄在0~10A范圍與電流關(guān)系滿足正弦函數(shù)，分別用0~5階樣條插值方法給出經(jīng)過(guò)數(shù)據(jù)點(diǎn)的數(shù)值逼近函數(shù)曲線。
Demo代碼
import matplotlib
import numpy as np
from matplotlib import pyplot as plt
from scipy import interpolate
font = {
"family": "Microsoft YaHei"
}
matplotlib.rc("font", **font)
# 創(chuàng)建數(shù)據(jù)點(diǎn)集
x = np.linspace(0, 10, 11)
y = np.sin(x)
# 繪制數(shù)據(jù)點(diǎn)集
plt.figure(figsize=(12, 9))
plt.plot(x, y, 'ro')
# 根據(jù)kind創(chuàng)建interp1d對(duì)象f、計(jì)算插值結(jié)果
xnew = np.linspace(0, 10, 101)
# 鄰接 0階 線性 二階
for kind in ['nearest', 'zero', 'linear', 'quadratic']:
f = interpolate.interp1d(x, y, kind=kind)
ynew = f(xnew)
plt.plot(xnew, ynew, label=str(kind))
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.legend(loc="lower right")
plt.show()
輸出：
分別對(duì)每一種插值方式進(jìn)行查看
1.當(dāng)kind = nearest時(shí)
import matplotlib
import numpy as np
from matplotlib import pyplot as plt
from scipy import interpolate
font = {
"family": "Microsoft YaHei"
}
matplotlib.rc("font", **font)
# 創(chuàng)建數(shù)據(jù)點(diǎn)集
x = np.linspace(0, 10, 11)
y = np.sin(x)
# 得插值函數(shù)
f = interpolate.interp1d(x, y, kind='nearest')
# 新數(shù)據(jù)
x_new = np.linspace(0,10,101)
y_new = f(x_new)
# 可視化
plt.plot(x, y, 'o', x_new, y_new, '-')
plt.show()
2.當(dāng)kind = zero時(shí)
import matplotlib
import numpy as np
from matplotlib import pyplot as plt
from scipy import interpolate
font = {
"family": "Microsoft YaHei"
}
matplotlib.rc("font", **font)
# 創(chuàng)建數(shù)據(jù)點(diǎn)集
x = np.linspace(0, 10, 11)
y = np.sin(x)
# 得插值函數(shù)
f = interpolate.interp1d(x, y, kind='zero')
# 新數(shù)據(jù)
x_new = np.linspace(0,10,101)
y_new = f(x_new)
# 可視化
plt.plot(x, y, 'o', x_new, y_new, '-')
plt.show()
3.當(dāng)kind = linear時(shí)
import matplotlib
import numpy as np
from matplotlib import pyplot as plt
from scipy import interpolate
font = {
"family": "Microsoft YaHei"
}
matplotlib.rc("font", **font)
# 創(chuàng)建數(shù)據(jù)點(diǎn)集
x = np.linspace(0, 10, 11)
y = np.sin(x)
# 得插值函數(shù)
f = interpolate.interp1d(x, y, kind='linear')
# 新數(shù)據(jù)
x_new = np.linspace(0,10,101)
y_new = f(x_new)
# 可視化
plt.plot(x, y, 'o', x_new, y_new, '-')
plt.show()
4.當(dāng)kind = quadratic時(shí)
import matplotlib
import numpy as np
from matplotlib import pyplot as plt
from scipy import interpolate
font = {
"family": "Microsoft YaHei"
}
matplotlib.rc("font", **font)
# 創(chuàng)建數(shù)據(jù)點(diǎn)集
x = np.linspace(0, 10, 11)
y = np.sin(x)
# 得插值函數(shù)
f = interpolate.interp1d(x, y, kind='quadratic')
# 新數(shù)據(jù)
x_new = np.linspace(0,10,101)
y_new = f(x_new)
# 可視化
plt.plot(x, y, 'o', x_new, y_new, '-')
plt.show()
5.當(dāng)kind = cubic時(shí)
import matplotlib
import numpy as np
from matplotlib import pyplot as plt
from scipy import interpolate
font = {
"family": "Microsoft YaHei"
}
matplotlib.rc("font", **font)
# 創(chuàng)建數(shù)據(jù)點(diǎn)集
x = np.linspace(0, 10, 11)
y = np.sin(x)
# 得插值函數(shù)
f = interpolate.interp1d(x, y, kind='cubic')
# 新數(shù)據(jù)
x_new = np.linspace(0,10,101)
y_new = f(x_new)
# 可視化
plt.plot(x, y, 'o', x_new, y_new, '-')
plt.show()
2. 二維插值
interp2d(x, y, z, kind=“'') 計(jì)算二維插值
2.1 圖像模糊處理——樣條插值
例3：某圖像表達(dá)式為，完成圖像的二維插值使其變清晰。
Demo代碼
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
def func(x, y):
return (x + y) * np.exp(-5.0 * (x ** 2 + y ** 2))
# X-Y軸分為15*15的網(wǎng)格
# x, y = np.mgrid[-1:1:15j, -1:1:15j]
x = np.linspace(-1, 1, 15)
y = np.linspace(-1, 1, 15)
x, y = np.meshgrid(x, y)
fvals = func(x, y)
# 二維插值
newfunc = interpolate.interp2d(x, y, fvals, kind='cubic')
# 計(jì)算100*100網(wǎng)格上插值
xnew = np.linspace(-1, 1, 100)
ynew = np.linspace(-1, 1, 100)
fnew = newfunc(xnew, ynew)
xnew, ynew = np.meshgrid(xnew, ynew)
plt.subplot(121)
# extent x軸和y軸范圍
im1 = plt.imshow(fvals, extent=[-1, 1, -1, 1], interpolation="nearest", origin="lower",cmap="Reds")
plt.colorbar(im1)
plt.subplot(122)
im2 = plt.imshow(fnew, extent=[-1, 1, -1, 1], interpolation="nearest", origin="lower",cmap="Reds")
plt.colorbar(im2)
plt.show()
輸出：
2.2 二維插值的三維圖
例4：某圖像表達(dá)式為，完成三維圖像的二維插值可視化。
?其實(shí)就是在二維插值基礎(chǔ)上 實(shí)現(xiàn)了三維圖像的繪制
Demo代碼
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
from scipy import interpolate
import matplotlib.cm as cm
import matplotlib.pyplot as plt
def func(x, y):
return (x + y) * np.exp(-5.0 * (x ** 2 + y ** 2))
# X-Y軸分為20*20的網(wǎng)格
x = np.linspace(-1, 1, 20)
y = np.linspace(-1, 1, 20)
x, y = np.meshgrid(x, y)
fvals = func(x, y)
# 繪制分圖1
fig = plt.figure(figsize=(9, 6))
ax = plt.subplot(1, 2, 1, projection='3d')
surf = ax.plot_surface(x, y, fvals, rstride=2, cstride=2, cmap=cm.coolwarm, linewidth=0.5, antialiased=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x,y)')
plt.colorbar(surf, shrink=0.5, aspect=5) # 添加顏色條標(biāo)注
# 二維插值
newfunc = interpolate.interp2d(x, y, fvals, kind='cubic')
# 計(jì)算100*100網(wǎng)格上插值
xnew = np.linspace(-1, 1, 100)
ynew = np.linspace(-1, 1, 100)
fnew = newfunc(xnew, ynew)
xnew, ynew = np.meshgrid(xnew, ynew)
ax2 = plt.subplot(1, 2, 2, projection='3d')
surf2 = ax2.plot_surface(xnew, ynew, fnew, rstride=2, cstride=2, cmap=cm.coolwarm, linewidth=0.5, antialiased=True)
ax2.set_xlabel('xnew')
ax2.set_ylabel('ynew')
ax2.set_zlabel('fnew(x,y)')
plt.colorbar(surf2, shrink=0.5, aspect=5)
# 標(biāo)注
plt.show()
輸出：
3. 最小二乘擬合
擬合指的是已知某函數(shù)的若干離散函數(shù)值{f1,f2,…,fn}，通過(guò)調(diào)整該函 數(shù)中若干待定系數(shù)f(λ1, λ2,…,λn)，使得該函數(shù)與已知點(diǎn)集的差別(最小二乘意義)最小。
如果待定函數(shù)是線性，就叫線性擬合或者線性回歸(主要在統(tǒng)計(jì)中)，否則叫作非線性擬合或者非線性回歸。表達(dá)式也可以是分段函數(shù)，這種情況下叫作樣條擬合。
從幾何意義上講，擬合是給定了空間中的一些點(diǎn)，找到一個(gè)已知形式、未知參數(shù)的連續(xù)曲面來(lái)最大限度地逼近這些點(diǎn)；而插值是找到一個(gè)（或幾個(gè)分片光滑的）連續(xù)曲面來(lái)穿過(guò)這些點(diǎn)。
選擇參數(shù)c使得擬合模型與實(shí)際觀測(cè)值在曲線擬合各點(diǎn)的殘差（或離差）ek=yk-f(xk,c)的加權(quán)平方和達(dá)到最小,此時(shí)所求曲線稱作在加權(quán)最小二乘意義下對(duì)數(shù)據(jù)的擬合曲線，這種方法叫做最小二乘法。
涉及知識(shí)點(diǎn)
from scipy.optimize import leastsq
例5：對(duì)下列電學(xué)元件的電壓電流記錄結(jié)果進(jìn)行最小二乘擬合，繪制相應(yīng)曲線。電流（A）8.19 2.72 6.39 8.71 4.7 2.66 3.78 電壓（V）7.01 2.78 6.47 6.71 4.1 4.23 4.05
Demo代碼
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
# 引入中文字體
font = {
"family": "Microsoft YaHei"
}
matplotlib.rc("font", **font)
# 設(shè)置圖字號(hào)
plt.figure(figsize=(9, 9))
# 初始數(shù)據(jù)值
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])
# 計(jì)算以p為參數(shù)的直線與原始數(shù)據(jù)之間誤差
def f(p):
k, b = p
return (Y - (k * X + b))
# leastsq使得f的輸出數(shù)組的平方和最小，參數(shù)初始值k、b設(shè)為[1,0]
r = leastsq(f, [1, 0])
# 得到計(jì)算出的最優(yōu)k、b
k, b = r[0]
# 可視化
plt.scatter(X, Y, s=100, alpha=1.0, marker='o', label=u'數(shù)據(jù)點(diǎn)')
x = np.linspace(0, 10, 1000)
y = k * x + b
ax = plt.gca()
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()
輸出：
結(jié)語(yǔ)
學(xué)習(xí)來(lái)源：B站及其課堂PPT，對(duì)其中代碼進(jìn)行了復(fù)現(xiàn)
鏈接：
https://www.bilibili.com/video/BV12h411d7Dm?from=search&seid=5685064698782810720
文章僅作為學(xué)習(xí)筆記，記錄從0到1的一個(gè)過(guò)程
希望對(duì)您有所幫助，如有錯(cuò)誤歡迎小伙伴指正～
評(píng)論
圖片
表情