1. 數(shù)據(jù)插值

插值是一種從離散數(shù)據(jù)點(diǎn)構(gòu)建函數(shù)的數(shù)學(xué)方法。插值函數(shù)或者插值方法應(yīng)該與給定的數(shù)據(jù)點(diǎn)完全一致。插值可能的應(yīng)用場景：

1. 根據(jù)給定的數(shù)據(jù)集繪制平滑的曲線
2. 對(duì)計(jì)算量很大的復(fù)雜函數(shù)進(jìn)行近似求值

插值和前面介紹過的最小二乘擬合有些類似。在最小二乘擬合中，我們感興趣的是使用數(shù)據(jù)點(diǎn)和超定方程組，將函數(shù)擬合到數(shù)組點(diǎn)，使得誤差平方和最小。在插值中，我們需要一個(gè)方程能夠與已有的數(shù)據(jù)點(diǎn)完全重合，僅使用與插值函數(shù)自由參數(shù)個(gè)數(shù)相同的數(shù)據(jù)點(diǎn)。因此，最小二乘法適合將大量數(shù)據(jù)點(diǎn)擬合到模型函數(shù)，插值是根據(jù)少量數(shù)據(jù)點(diǎn)創(chuàng)建函數(shù)。

外插（extrapolation）是與插值（interpolation）相關(guān)的一個(gè)概念。外插是在采樣范圍之外計(jì)算函數(shù)的估計(jì)值。我們這里只介紹插值。

2. 導(dǎo)入模塊

本部分我們將使用NumPy中的polynomial模塊和SciPy的interpolation模塊。

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt

from numpy import polynomial as P
from scipy import interpolate
%reload_ext version_information
%version_information numpy, matplotlib, scipy

3. 插值函數(shù)

為了簡潔起見，我們這里只考慮一維插值問題。對(duì)于給定的數(shù)據(jù)點(diǎn)的集合 { ( x i , y i ) } i = 1 n \left\{ (x_i, y_i) \right \}_{i=1}^n {(xi?,yi?)}i=1n?，找到插值函數(shù)$f(x_i)=y_i$。插值函數(shù)的選擇并不是唯一的，事實(shí)上有無數(shù)函數(shù)滿足插值標(biāo)準(zhǔn)。通?？梢园巡逯岛瘮?shù)寫為一組基函數(shù)$?(x)$的線性組合，即$f(x)=\sum _{j=1}^n{c}_j{?}_j(x)$，其中$c_j$是未知系數(shù)。將給定的數(shù)據(jù)點(diǎn)代入線性組合，可以得到未知系數(shù)的線性方程組：

$$?\begin{array}{cccc}{?}_1(x_1)& {?}_2(x_1)& ?& {?}_n(x_1)\\ {?}_1(x_2)& {?}_2(x_2)& ?& {?}_n(x_2)\\ ?& ?& ?& ?\\ {?}_1(x_n)& {?}_2(x_n)& ?& {?}_n(x_n)\end{array}??\begin{array}{c}{c}_1\\ c_2\\ ?\\ c_n\end{array}??\begin{array}{c}{y}_1\\ y_2\\ ?\\ y_n\end{array}?$$

這里基函數(shù)的數(shù)量與數(shù)據(jù)點(diǎn)的數(shù)量想用，可以使用求解方程組的標(biāo)準(zhǔn)方法得到系數(shù)向$c$的唯一解。如果基函數(shù)的數(shù)量小于數(shù)據(jù)點(diǎn)的數(shù)量，該方程組是超定的，可能需要使用最小二乘擬合進(jìn)行近似插值。

3.1 多項(xiàng)式

NumPy庫包含的polynominal模塊提供了處理多項(xiàng)式的函數(shù)和類。要?jiǎng)?chuàng)建Polynominal類的實(shí)例，可以將系數(shù)數(shù)組傳給該類的構(gòu)造函數(shù)。

p1 = P.Polynomial([1,2,3])
p1

也可以使用Polynomial.fromroots方法，通過指定多項(xiàng)式的根來初始化多項(xiàng)式。

p2 = P.Polynomial.fromroots([-1, 1])
p2

我們可以通過實(shí)例的方法訪問其特性：

p1.roots() # 根

p1.coef # 系數(shù)

Polynomial實(shí)例包括domain和windows兩個(gè)參數(shù)，可以用于把多項(xiàng)式的輸入域映射到另外一個(gè)區(qū)間。

p1.domain

p1.window

多項(xiàng)式在用Polynomial實(shí)例表示后可以很容易計(jì)算任何$x$的多項(xiàng)式值。

p1(np.array([1.5, 2.5, 3.5]))

Polynomial實(shí)例支持標(biāo)準(zhǔn)的算術(shù)運(yùn)算符。

p1 + p2

分解因式

p3 = P.Polynomial([1,1])  # 等價(jià)于 P.Polynomial.fromroots(-1)
p2 // p3

P.Polynomial.fromroots(-1)

除了標(biāo)準(zhǔn)冪基多項(xiàng)式的Polynomial類，polynomial模塊還提供了切比雪夫多項(xiàng)式、勒讓德多項(xiàng)式、拉蓋爾多項(xiàng)式、埃爾米特多項(xiàng)式的類。

l1 = P.Legendre([1,2,3])  # 勒讓德多項(xiàng)式
l1

l1.roots()

3.2 多項(xiàng)式插值

求解多項(xiàng)式插值問題，首先需要根據(jù)提供的基函數(shù)計(jì)算矩陣$?(x)$，然后求解得到的線性方程組。polynomial模塊中的每個(gè)多項(xiàng)式類都提供了方便的函數(shù)來計(jì)算相應(yīng)基的范德蒙矩陣。

例如，假設(shè)存在數(shù)據(jù)點(diǎn)(1, 1)、(2, 3)、(3, 5)、(4, 4)。

x = np.array([1, 2, 3, 4])
y = np.array([1, 3, 5, 4])
deg = len(x) - 1
A = P.polynomial.polyvander(x, deg)  # 范德蒙矩陣
A

c = linalg.solve(A, y)
c

找到的插值多項(xiàng)式是$f(x)=2?3.5x+3x^2?0.5x^3$。

f1 = P.Polynomial(c)
f1(2.5)

下面，我們使用切比雪夫多項(xiàng)式重新對(duì)數(shù)據(jù)進(jìn)行插值。

A = P.chebyshev.chebvander(x, deg)
A

c = linalg.solve(A, y)
c

f2 = P.Chebyshev(c)
f2(2.5)

將函數(shù)$f1$和$f2$進(jìn)行對(duì)比，驗(yàn)證不同基函數(shù)進(jìn)行插值，得到了一致的插值函數(shù)。

xx = np.linspace(x.min(), x.max(), 100)

fig, ax = plt.subplots(1, 1, figsize=(8, 4))

ax.plot(xx, f1(xx), 'b', lw=2, label='Power basis interp.')
ax.plot(xx, f2(xx), 'r--', lw=2, label='Chebyshev basis interp.')
ax.scatter(x, y, label='data points')

ax.legend(loc=4)
ax.set_xticks(x)
ax.set_ylabel(r"$y$", fontsize=18)
ax.set_xlabel(r"$x$", fontsize=18)

polynomial模塊還提供了更快捷的方法計(jì)算插值多項(xiàng)式。每個(gè)多項(xiàng)式類都有fit方法用于計(jì)算插值，deg參數(shù)用于控制多項(xiàng)式的次數(shù)。如果多項(xiàng)式的次數(shù)少于數(shù)據(jù)點(diǎn)的數(shù)目減1，fit方法會(huì)使用最小二乘擬合而不是精確插值。

f1b = P.Polynomial.fit(x, y, deg)
f1b

f2b = P.Chebyshev.fit(x, y, deg)
f2b

注意在使用fit方法時(shí)，實(shí)例中的domain屬性會(huì)自動(dòng)設(shè)置為數(shù)據(jù)點(diǎn)中相應(yīng)的x值，上述示例中為[1, 4]，并相應(yīng)調(diào)整系數(shù)。將插值數(shù)據(jù)映射到某個(gè)基函數(shù)最合適的范圍，可以明顯提高插值的數(shù)值穩(wěn)定性，減小范德蒙矩陣的條件數(shù)。

np.linalg.cond(P.chebyshev.chebvander(x, deg))

np.linalg.cond(P.chebyshev.chebvander((2*x-5)/3.0, deg))

當(dāng)數(shù)據(jù)點(diǎn)的數(shù)量增加時(shí)，需要使用次數(shù)更多的多項(xiàng)式才能得到精確的插值。這會(huì)帶來多方面問題：

1. 次數(shù)更多的多項(xiàng)式插值計(jì)算和插值函數(shù)的確定，計(jì)算量很大。
2. 高次多項(xiàng)式插值可能會(huì)在插值點(diǎn)之間帶來不可預(yù)料的行為，插值函數(shù)在數(shù)據(jù)點(diǎn)之間可能變化很大。
   下面我們將使用龍格函數(shù)演示這種行為：

def runge(x):
    return 1/(1 + 25 * x**2)

def runge_interpolate(n):
    x = np.linspace(-1, 1, n+1)
    p = P.Polynomial.fit(x, runge(x), deg=n)
    return x, p

xx = np.linspace(-1, 1, 250)

fig, ax = plt.subplots(1, 1, figsize=(8, 4))

ax.plot(xx, runge(xx), 'k', lw=2, label="Runge's function")

n = 13
x, p = runge_interpolate(n)
ax.plot(x, runge(x), 'ro')
ax.plot(xx, p(xx), 'r', label='interp. order %d' % n)

n = 14
x, p = runge_interpolate(n)
ax.plot(x, runge(x), 'go')
ax.plot(xx, p(xx), 'g', label='interp. order %d' % n)

ax.legend(loc=8)
ax.set_xlim(-1.1, 1.1)
ax.set_ylim(-1, 2)
ax.set_xticks([-1, -0.5, 0, 0.5, 1])
ax.set_ylabel(r"$y$", fontsize=18)
ax.set_xlabel(r"$x$", fontsize=18)

可以看到，高階插值函數(shù)在遠(yuǎn)端采樣點(diǎn)之間急劇震蕩。這種不良特性違背了插值的初衷。一個(gè)可行的解決方法是，當(dāng)面對(duì)大量數(shù)據(jù)點(diǎn)時(shí)，使用分段低次多項(xiàng)式進(jìn)行插值。

3.3 樣條插值

對(duì)于包含$n$個(gè)數(shù)據(jù)點(diǎn)的集合，在整個(gè)數(shù)據(jù)區(qū)間上有$n?1$個(gè)子區(qū)間。連接兩個(gè)子區(qū)間的數(shù)據(jù)點(diǎn)在分段多項(xiàng)式插值中被稱為節(jié)點(diǎn)。在每個(gè)子區(qū)間上使用$k$次多項(xiàng)式對(duì)$n$個(gè)數(shù)據(jù)點(diǎn)進(jìn)行插值，需要確定$(k+1)(n?1)$個(gè)參數(shù)。所有節(jié)點(diǎn)的值可以給出$2(n?1)$個(gè)方程。為了保證分段多項(xiàng)式的平滑，節(jié)點(diǎn)處導(dǎo)數(shù)和高階導(dǎo)數(shù)的連續(xù)性也會(huì)給出相應(yīng)方程。

樣條是一種特殊類型的分段式插值函數(shù)。最常用的是三次樣條，$k=3$，需要$4(n?1)$個(gè)參數(shù)。在$n?2$個(gè)節(jié)點(diǎn)處，一階和二階導(dǎo)數(shù)的連續(xù)性可以給出$2(n?1)$個(gè)方程，總方程的數(shù)目為$4(n?1)?2$。此時(shí)還剩下兩個(gè)未確定的參數(shù)。一種常見的方法是要求端點(diǎn)處的二階導(dǎo)數(shù)為0，這被稱為自然樣條。

SciPy的interpolate模塊提供了用于樣條插值的多個(gè)函數(shù)和類。下面我們將再次使用龍格函數(shù)，演示使用interpolate.interp1d函數(shù)。這里設(shè)置kind=3來計(jì)算三次樣條。

x = np.linspace(-1, 1, 11)
y = runge(x)
f = interpolate.interp1d(x, y, kind=3)
f(0.05)

xx = np.linspace(-1, 1, 100)
fig, ax = plt.subplots(figsize=(8, 4))

ax.plot(xx, runge(xx), 'k', lw=1, label="Runge's function")
ax.plot(x, y, 'ro', label='sample points')
ax.plot(xx, f(xx), 'r--', lw=2, label='spline order 3')

ax.legend()
ax.set_ylim(0, 1.1)
ax.set_xticks([-1, -0.5, 0, 0.5, 1])
ax.set_ylabel(r"$y$", fontsize=18)
ax.set_xlabel(r"$x$", fontsize=18)

這里使用了11個(gè)數(shù)據(jù)點(diǎn)和三次樣條，可以看到插值很好地與原始函數(shù)吻合。
為了說明階數(shù)對(duì)樣條插值的影響，我們下面準(zhǔn)備了8個(gè)數(shù)據(jù)點(diǎn)使用不同階數(shù)的樣條進(jìn)行插值。

x = np.array([0, 1, 2, 3, 4, 5, 6, 7])
y = np.array([3, 4, 3.5, 2, 1, 1.5, 1.25, 0.9])

xx = np.linspace(x.min(), x.max(), 100)

fig, ax = plt.subplots(figsize=(8, 4))

ax.scatter(x, y)

for n in [1, 2, 3, 5]:
    f = interpolate.interp1d(x, y, kind=n)
    ax.plot(xx, f(xx), label='order %d' % n)

ax.legend()
ax.set_ylabel(r"$y$", fontsize=18)
ax.set_xlabel(r"$x$", fontsize=18)

可以看到，二階或三階樣條已經(jīng)能夠提供較好的插值。高階樣條會(huì)出現(xiàn)高階多項(xiàng)式插值中類似的數(shù)值震蕩問題。

3.4 多變量插值

多項(xiàng)式插值和樣條插值都可以直接推廣到多變量情況。SciPy為多變量插值提供了多個(gè)函數(shù)和類。我們將介紹兩個(gè)最常用的雙變量插值函數(shù)：interpolate.interp2d和interpolate.griddata。

3.4.1 均勻網(wǎng)格

interpolate.interp2d函數(shù)是interpolate.interp1d函數(shù)的推廣，要求數(shù)據(jù)點(diǎn)位于x和y坐標(biāo)組成的規(guī)則且均勻的網(wǎng)格中。

為了演示函數(shù)用法，我們?cè)谝阎瘮?shù)
$$f(x,y)=\exp (?(x+1/2{)}^2?2(y+1/2{)}^2)?\exp (?(x?1/2{)}^2?2(y?1/2{)}^2)$$
中添加隨機(jī)噪聲來模擬測量噪聲。為了構(gòu)造插值問題，我們?cè)趚和y坐標(biāo)的[-2, 2]區(qū)間上采樣10個(gè)點(diǎn)。

def f(x, y):
    return np.exp(-(x + .5)**2 - 2*(y + .5)**2) - np.exp(-(x - .5)**2 - 2*(y - .5)**2)

x = y = np.linspace(-2, 2, 10)
X, Y = np.meshgrid(x, y)
# simulate noisy data at fixed grid points X, Y
Z = f(X, Y) + 0.05  * np.random.randn(*X.shape)
f_interp = interpolate.interp2d(x, y, Z, kind='cubic')  # 三次樣條插值
xx = yy = np.linspace(x.min(), x.max(), 100)
ZZ_interp = f_interp(xx, yy)
XX, YY = np.meshgrid(xx, yy)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

c = axes[0].contourf(XX, YY, f(XX, YY), 15, cmap=plt.cm.RdBu)
axes[0].set_xlabel(r"$x$", fontsize=20)
axes[0].set_ylabel(r"$y$", fontsize=20)
axes[0].set_title("exact / high sampling")
cb = fig.colorbar(c, ax=axes[0])
cb.set_label(r"$z$", fontsize=20)

c = axes[1].contourf(XX, YY, ZZ_interp, 15, cmap=plt.cm.RdBu)
axes[1].set_ylim(-2.1, 2.1)
axes[1].set_xlim(-2.1, 2.1)
axes[1].set_xlabel(r"$x$", fontsize=20)
axes[1].set_ylabel(r"$y$", fontsize=20)
axes[1].scatter(X, Y, marker='x', color='k')
axes[1].set_title("interpolation of noisy data / low sampling")
cb = fig.colorbar(c, ax=axes[1])
cb.set_label(r"$z$", fontsize=20)

對(duì)于更高維的問題，可以使用interpolate.interpnd函數(shù)，它是對(duì)N維問題的推廣。

3.4.2 不均勻網(wǎng)格

另一種多變量插值的常見場景是從不規(guī)則的坐標(biāo)網(wǎng)絡(luò)中采樣數(shù)據(jù)，例如實(shí)驗(yàn)室數(shù)據(jù)等。為了能夠使用現(xiàn)有的工具對(duì)數(shù)據(jù)進(jìn)行可視化分析，可能需要將它們插值到規(guī)則的坐標(biāo)網(wǎng)格中。

在SciPy中，可以使用interpolate.griddata函數(shù)來完成這項(xiàng)工作。該函數(shù)的第一個(gè)參數(shù)是坐標(biāo)向量構(gòu)成的元組，第二個(gè)參數(shù)是對(duì)應(yīng)的數(shù)據(jù)點(diǎn)的值，第三個(gè)參數(shù)是待插值的新數(shù)據(jù)點(diǎn)的坐標(biāo)數(shù)組或坐標(biāo)矩陣。
我們考慮函數(shù)$f(x,y)=\exp (?x^2?y^2)\cos 4x\sin 6y$在x和y坐標(biāo)區(qū)間[-1, 1]上隨機(jī)選擇的采樣點(diǎn)，并對(duì)采樣數(shù)據(jù)進(jìn)行插值。

def f(x, y):
    return np.exp(-x**2 - y**2) * np.cos(4*x) * np.sin(6*y)

x = y = np.linspace(-1, 1, 100)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)
np.random.seed(0)
N = 500
xdata = np.random.uniform(-1, 1, N)
ydata = np.random.uniform(-1, 1, N)
zdata = f(xdata, ydata)

我們畫出函數(shù) ??(??,??) 的等高線圖，以及500個(gè)采樣點(diǎn)的位置。

fig, ax = plt.subplots(figsize=(8, 6))
c = ax.contourf(X, Y, Z, 15, cmap=plt.cm.RdBu);
ax.scatter(xdata, ydata, marker='.')
ax.set_ylim(-1,1)
ax.set_xlim(-1,1)
ax.set_xlabel(r"$x$", fontsize=20)
ax.set_ylabel(r"$y$", fontsize=20)

cb = fig.colorbar(c, ax=ax)
cb.set_label(r"$z$", fontsize=20)

隨后，我們使用X和Y坐標(biāo)數(shù)組定義更細(xì)小間隔（超采樣）的網(wǎng)格中對(duì)數(shù)據(jù)進(jìn)行插值。我們將比較最近鄰點(diǎn)插值（0階插值）、線性插值（1階插值）和三次樣條插值在不用數(shù)據(jù)點(diǎn)數(shù)目下的效果區(qū)別。

def z_interpolate(xdata, ydata, zdata):
    Zi_0 = interpolate.griddata((xdata, ydata), zdata, (X, Y), method='nearest')
    Zi_1 = interpolate.griddata((xdata, ydata), zdata, (X, Y), method='linear')
    Zi_3 = interpolate.griddata((xdata, ydata), zdata, (X, Y), method='cubic')
    return Zi_0, Zi_1, Zi_3

fig, axes = plt.subplots(3, 3, figsize=(12, 12), sharex=True, sharey=True)

n_vec = [50, 150, 500]

for idx, n in enumerate(n_vec):
    Zi_0, Zi_1, Zi_3 = z_interpolate(xdata[:n], ydata[:n], zdata[:n])
    axes[idx, 0].contourf(X, Y, Zi_0, 15, cmap=plt.cm.RdBu)
    axes[idx, 0].set_ylabel("%d data points\ny" % n, fontsize=16)
    axes[idx, 0].set_title("nearest", fontsize=16)
    axes[idx, 1].contourf(X, Y, Zi_1, 15, cmap=plt.cm.RdBu)
    axes[idx, 1].set_title("linear", fontsize=16)
    axes[idx, 2].contourf(X, Y, Zi_3, 15, cmap=plt.cm.RdBu)
    axes[idx, 2].set_title("cubic", fontsize=16)

for m in range(len(n_vec)):
    axes[idx, m].set_xlabel("x", fontsize=16)

只要采樣點(diǎn)能夠有效覆蓋我們感興趣的區(qū)域，就可以通過非結(jié)構(gòu)化的樣本數(shù)據(jù)進(jìn)行插值來重建函數(shù)。上面例子中，三次樣條插值的效果要遠(yuǎn)優(yōu)于最近鄰點(diǎn)插值和線性插值。

參考文獻(xiàn)來自桑鴻乾老師的課件：科學(xué)計(jì)算和人工智能文章來源地址http://www.zghlxwxcb.cn/news/detail-598999.html

到了這里，關(guān)于【Python數(shù)據(jù)插值】的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！