平滑樣條

1D 中的樣條平滑

對於內插問題，任務是建構一條通過給定資料點集的曲線。如果資料有雜訊，這可能不適合：我們那時希望建構一條平滑曲線 \(g(x)\)，它近似輸入資料，而無需完全通過每個點。

為此，scipy.interpolate 允許建構平滑樣條，以平衡結果曲線 \(g(x)\) 與資料的接近程度，以及 \(g(x)\) 的平滑度。在數學上，任務是解決懲罰最小平方問題，其中懲罰控制 \(g(x)\) 的平滑度。

我們提供了兩種建構平滑樣條的方法，它們在 (1) 懲罰項的形式和 (2) 建構平滑曲線的基底方面有所不同。下面我們考慮這兩種方法。

前一種變體由 make_smoothing_spline 函數執行，它是 H. Woltring 經典 gcvspl Fortran 套件的潔淨室重新實作。後一種變體由 make_splrep 函數實作，它是 P. Dierckx 的 Fortran FITPACK 函式庫的重新實作。FITPACK 函式庫的舊版介面 也可用。

「經典」平滑樣條和廣義交叉驗證 (GCV) 準則

給定資料陣列 x 和 y 以及非負權重陣列 w，我們尋找三次樣條函數 g(x)，使其最小化

\[\sum\limits_{j=1}^n w_j \left\lvert y_j - g(x_j) \right\rvert^2 + \lambda\int\limits_{x_1}^{x_n} \left( g^{(2)}(u) \right)^2 d u\]

其中 \(\lambda \geqslant 0\) 是一個非負懲罰參數，而 \(g^{(2)}(x)\) 是 \(g(x)\) 的二階導數。第一項中的總和遍歷資料點 \((x_j, y_j)\)，而第二項中的積分遍歷整個區間 \(x \in [x_1, x_n]\)。

這裡，第一項懲罰樣條函數與資料的偏差，第二項懲罰二階導數的大值—這被視為曲線平滑度的準則。

目標函數 \(g(x)\) 被視為以資料點 \(x_j\) 為節點的自然三次樣條，並且在給定 \(\lambda\) 值下對樣條係數進行最小化。

顯然，\(\lambda = 0\) 對應於內插問題—結果是自然內插樣條；在相反的極限下，\(\lambda \gg 1\)，結果 \(g(x)\) 接近一條直線（因為最小化有效地將 \(g(x)\) 的二階導數歸零）。

平滑函數強烈依賴於 \(\lambda\)，並且有多種策略可用於選擇「最佳」懲罰值。一種流行的策略是所謂的廣義交叉驗證 (GCV)：從概念上講，這相當於比較在簡化資料集上建構的樣條函數，其中我們省略了一個資料點。直接應用這種留一法交叉驗證程序非常耗費資源，我們使用更有效率的 GCV 演算法。

為了建構給定資料和懲罰參數的平滑樣條，我們使用函數 make_smoothing_spline。它的介面類似於內插樣條的建構子 make_interp_spline：它接受資料陣列並傳回可呼叫的 BSpline 實例。

此外，它接受一個可選的 lam 關鍵字引數來指定懲罰參數 \(\lambda\)。如果省略或設定為 None，則 \(\lambda\) 會透過 GCV 程序計算。

為了說明懲罰參數的效果，請考慮一個帶有一些雜訊的正弦曲線的玩具範例

>>> import numpy as np
>>> from scipy.interpolate import make_smoothing_spline

產生一些雜訊資料

>>> x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/16)
>>> rng = np.random.default_rng()
>>> y =  np.sin(x) + 0.4*rng.standard_normal(size=len(x))

建構並繪製一系列懲罰參數值的平滑樣條

>>> import matplotlib.pyplot as plt
>>> xnew = np.arange(0, 9/4, 1/50) * np.pi
>>> for lam in [0, 0.02, 10, None]:
...     spl = make_smoothing_spline(x, y, lam=lam)
...     plt.plot(xnew, spl(xnew), '-.', label=fr'$\lambda=${lam}')
>>> plt.plot(x, y, 'o')
>>> plt.legend()
>>> plt.show()

我們清楚地看到 lam=0 建構了內插樣條；lam 的大值將結果曲線展平為一條直線；而 GCV 結果 lam=None 接近底層的正弦曲線。

具有自動節點選擇的平滑樣條

作為 make_smoothing_spline 的補充，SciPy 提供了另一種替代方案，形式為 make_splrep 和 make_splprep 例程。前者建構樣條函數，後者用於 \(d > 1\) 維度的參數樣條曲線。

雖然具有類似的 API（接收資料陣列，傳回 BSpline 實例），但這些與 make_smoothing_spline 在以下幾個方面有所不同

• 懲罰項的函數形式不同：這些例程使用 \(k\) 階導數的跳躍，而不是 \((k-1)\) 階導數的積分；

• 使用平滑參數 \(s\)，而不是懲罰參數 \(\lambda\)；

• 這些例程自動建構節點向量；根據輸入，結果樣條可能具有比資料點少得多的節點。

• 預設情況下，邊界條件不同：雖然 make_smoothing_spline 建構自然三次樣條，但這些例程預設使用非節點邊界條件。

讓我們更詳細地考慮該演算法。首先，是平滑準則。給定一個由節點 \(t_j\) 和係數 \(c_j\) 定義的三次樣條函數 \(g(x)\)，考慮內部節點處三階導數的跳躍，

\[D_j = g^{(3)}(t_j + 0) - g^{(3)}(t_j - 0) .\]

（對於 \(k\) 階樣條，我們將使用 \(k\) 階導數的跳躍。）

如果所有 \(D_j = 0\)，則 \(g(x)\) 是由節點跨越的整個域上的一個單一多項式。因此，我們認為 \(g(x)\) 是一個分段 \(C^2\) 可微樣條函數，並使用跳躍總和作為平滑準則，

\[\sum_j \left| D_j \right|^2 \to \mathrm{min} ,\]

其中最小化是在樣條係數以及可能的樣條節點上執行的。

為了確保 \(g(x)\) 近似輸入資料 \(x_j\) 和 \(y_j\)，我們引入平滑參數 \(s \geqslant 0\)，並新增約束條件，即

\[\sum_{j=1}^m \left[w_j \times \left(g(x_j) - y_j\right) \right]^2 \leqslant s .\]

在此公式中，平滑參數 \(s\) 是使用者輸入，很像經典平滑樣條的懲罰參數 \(\lambda\)。

請注意，極限 s = 0 對應於內插問題，其中 \(g(x_j) = y_j\)。增加 s 會導致更平滑的擬合，並且在非常大的 s 的極限下，\(g(x)\) 會退化為單一最佳擬合多項式。

對於固定的節點向量和給定的 \(s\) 值，最小化問題是線性的。如果我們也針對節點進行最小化，則問題會變成非線性。因此，我們需要指定一個迭代最小化過程，以建構節點向量以及樣條係數。

因此，我們使用以下程序

• 我們從沒有內部節點的樣條開始，並檢查使用者提供的 \(s\) 值的平滑條件。如果滿足，我們就完成了。否則，

• 迭代，其中在每次迭代中，我們

  • 透過分割樣條函數 \(g(x_j)\) 和資料 \(y_j\) 之間偏差最大的區間來新增新的節點。

  • 建構 \(g(x)\) 的下一個近似值，並檢查平滑準則。

如果滿足平滑條件，或達到允許的最大節點數，則迭代停止。後者可以由使用者指定，或作為預設值 len(x) + k + 1，這對應於使用 k 階樣條內插資料陣列 x。

重新措辭並忽略細節，該程序是在由 generate_knots 產生的節點向量上迭代，在每個步驟應用 make_lsq_spline。在虛擬碼中

for t in generate_knots(x, y, s=s):
    g = make_lsq_spline(x, y, t=t)     # construct
    if ((y - g(x))**2).sum() < s:      # check smoothness
        break

注意

對於 s=0，我們採用捷徑，並使用非節點邊界條件建構內插樣條，而不是迭代。

建構節點向量的迭代程序可透過產生器函數 generate_knots 取得。為了說明

>>> import numpy as np
>>> from scipy.interpolate import generate_knots
>>> x = np.arange(7)
>>> y = x**4
>>> list(generate_knots(x, y, s=1))   # default is cubic splines, k=3
[array([0., 0., 0., 0., 6., 6., 6., 6.]),
 array([0., 0., 0., 0., 3., 6., 6., 6., 6.]),
 array([0., 0., 0., 0., 3., 5., 6., 6., 6., 6.]),
 array([0., 0., 0., 0., 2., 3., 4., 6., 6., 6., 6.])]

對於 s=0，產生器會縮短

>>> list(generate_knots(x, y, s=0))
[array([0, 0, 0, 0, 2, 3, 4, 6, 6, 6, 6])]

注意

一般來說，節點放置在資料站點。偶數階樣條是例外，其中節點可以放置在遠離資料的位置。當 s=0（內插）時，或當 s 夠小時，會發生這種情況，以至於達到最大節點數，並且例程切換到 s=0 節點向量（有時稱為 Greville 橫座標）。

>>> list(generate_knots(x, y, s=1, k=2))   # k=2, quadratic spline
[array([0., 0., 0., 6., 6., 6.]),
 array([0., 0., 0., 3., 6., 6., 6.]),
 array([0., 0., 0., 3., 5., 6., 6., 6.]),
 array([0., 0., 0., 2., 3., 5., 6., 6., 6.]),
 array([0. , 0. , 0. , 1.5, 2.5, 3.5, 4.5, 6. , 6. , 6. ])]   # Greville sites

注意

建構節點向量的啟發式方法遵循 FITPACK Fortran 函式庫使用的演算法。該演算法是相同的，並且由於浮點捨入，可能會出現小的差異。

現在我們使用與上一節相同的玩具資料集來說明 make_splrep 結果

>>> import numpy as np
>>> from scipy.interpolate import make_splrep

產生一些雜訊資料

>>> x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/16)
>>> rng = np.random.default_rng()
>>> y =  np.sin(x) + 0.4*rng.standard_normal(size=len(x))

建構並繪製一系列 s 參數值的平滑樣條

>>> import matplotlib.pyplot as plt
>>> xnew = np.arange(0, 9/4, 1/50) * np.pi

>>> plt.plot(xnew, np.sin(xnew), '-.', label='sin(x)')
>>> plt.plot(xnew, make_splrep(x, y, s=0)(xnew), '-', label='s=0')
>>> plt.plot(xnew, make_splrep(x, y, s=len(x))(xnew), '-', label=f's={len(x)}')
>>> plt.plot(x, y, 'o')
>>> plt.legend()
>>> plt.show()

我們看到 \(s=0\) 曲線遵循資料點的（隨機）波動，而 \(s > 0\) 曲線接近底層的正弦函數。另請注意，外插值根據 \(s\) 的值而發生巨大變化。

找到 \(s\) 的良好值是一個反覆試驗的過程。如果權重對應於輸入資料標準差的倒數，則 \(s\) 的「良好」值預期會在 \(m - \sqrt{2m}\) 和 \(m + \sqrt{2m}\) 之間，其中 \(m\) 是資料點的數量。如果所有權重都等於 1，則合理的選擇可能在 \(s \sim m\,\sigma^2\) 附近，其中 \(\sigma\) 是資料標準差的估計值。

注意

節點數非常強烈地依賴於 s。小的 s 變化可能會導致節點數的劇烈變化。

注意

make_smoothing_spline 和 make_splrep 都允許加權擬合，其中使用者提供權重陣列 w。請注意，定義略有不同：make_smoothing_spline 將權重平方以與 gcvspl 保持一致，而 make_splrep 不會—以與 FITPACK 保持一致。

\(d>1\) 中的平滑樣條曲線

到目前為止，我們考慮了建構平滑樣條函數 \(g(x)\)，給定資料陣列 x 和 y。現在我們考慮建構平滑樣條曲線的相關問題，其中我們將資料視為平面上的點 \(\mathbf{p}_j = (x_j, y_j)\)，並且我們想要建構一個參數函數 \(\mathbf{g}(\mathbf{p}) = (g_x(u), g_y(u))\)，其中參數 \(u_j\) 的值對應於 \(x_j\) 和 \(y_j\)。

請注意，此問題很容易推廣到更高的維度 \(d > 2\)：我們只需要 \(d\) 個資料陣列，並建構一個具有 \(d\) 個分量的參數函數。

另請注意，參數化的選擇無法自動化，即使對於 內插曲線，不同的參數化也可能導致相同資料的曲線非常不同。

一旦選擇了特定的參數化形式，建構平滑曲線的問題在概念上與建構平滑樣條函數非常相似。簡而言之，

• 從參數 \(u\) 的值新增樣條節點，並且

• 我們為 樣條函數 考慮的要最小化的成本函數和約束條件只是在 \(d\) 個分量上新增了一個額外的總和。

make_splrep 函數的「參數」推廣是 make_splprep，其文件字串詳細說明了它解決的最小化問題的精確數學公式。

參數情況的主要使用者可見差異是用戶介面

• 不是兩個資料陣列 x 和 y，make_splprep 接收一個二維陣列，其中第二個維度的大小為 \(d\)，並且每個資料陣列都沿著第一個維度儲存（或者，您可以提供 1D 陣列的清單）。

• 傳回值是一對：一個 BSpline 實例和參數值陣列 u，它對應於輸入資料陣列。

預設情況下，make_splprep 建構並傳回輸入資料的弦長參數化（有關詳細資訊，請參閱 參數樣條曲線 區段）。或者，您可以提供您自己的參數值陣列 u。

為了說明 API，請考慮一個玩具問題：我們有一些從笛卡爾葉形線加上一些雜訊取樣的資料。

>>> import numpy as np
>>> from scipy.interpolate import make_splprep
>>> th = np.linspace(-0.2, np.pi/2 + 0.2, 21)
>>> r = 3 * np.sin(th) * np.cos(th) / (np.sin(th)**3 + np.cos(th)**3)
>>> x, y = r * np.cos(th), r * np.sin(th)

新增一些雜訊並建構內插器

>>> rng = np.random.default_rng()
>>> xn = x + 0.1*rng.uniform(-1, 1, size=len(x))
>>> yn = y + 0.1*rng.uniform(-1, 1, size=len(x))
>>> spl, u = make_splprep([xn, yn], s=0)   # note the [xn, yn] argument
>>> spl_n, u_n = make_splprep([xn, yn], s=0.1)

並繪製結果（spl(u) 的結果是一個 2D 陣列，因此我們將其解壓縮為一對 x 和 y 陣列以進行繪圖）。

>>> import matplotlib.pyplot as plt
>>> plt.plot(xn, yn, 'o')
>>> plt.plot(*spl(u), '--')
>>> plt.plot(*spl_n(u_n), '-')
>>> plt.show()

1-D 樣條平滑的舊版例程

除了我們在前面章節中討論的平滑樣條建構子之外，scipy.interpolate 還為 P. Dierckx 撰寫的歷史悠久的 FITPACK Fortran 函式庫提供了直接介面。

注意

這些介面應被視為舊版—雖然我們不打算棄用或移除它們，但我們建議新程式碼改用更現代的替代方案 make_smoothing_spline、make_splrep 或 make_splprep。

由於歷史原因，scipy.interpolate 為 FITPACK 提供了兩個等效的介面，一個程序性介面和一個物件導向介面。雖然等效，但這些介面具有不同的預設值。下面我們將依次討論它們，從功能性介面開始。

程序性 (splrep)

樣條內插需要兩個基本步驟：(1) 計算曲線的樣條表示，以及 (2) 在所需點評估樣條。為了找到樣條表示，有兩種不同的方法來表示曲線並獲得（平滑）樣條係數：直接和參數化。直接方法使用函數 splrep 找到 2-D 平面中曲線的樣條表示。前兩個引數是唯一需要的，它們提供了曲線的 \(x\) 和 \(y\) 分量。正常輸出是一個 3 元組 \(\left(t,c,k\right)\)，其中包含節點點 \(t\)、係數 \(c\) 和樣條的階數 \(k\)。預設樣條階數為三次，但可以使用輸入關鍵字 k 更改。

節點陣列將內插區間定義為 t[k:-k]，因此 t 陣列的前 \(k+1\) 個和最後 \(k+1\) 個條目定義了邊界節點。係數是一個 1D 陣列，長度至少為 len(t) - k - 1。某些例程會填充此陣列，使其具有 len(c) == len(t)—這些額外係數在樣條評估中會被忽略。

tck 元組格式與 內插 b 樣條 相容：splrep 的輸出可以包裝到 BSpline 物件中，例如 BSpline(*tck)，並且下面描述的評估/積分/求根例程可以互換使用 tck 元組和 BSpline 物件。

對於 N-D 空間中的曲線，函數 splprep 允許參數化定義曲線。對於此函數，只需要 1 個輸入引數。此輸入是表示 N-D 空間中曲線的 \(N\) 個陣列的清單。每個陣列的長度是曲線點的數量，每個陣列提供 N-D 資料點的一個分量。參數變數由關鍵字引數 u 給出，預設為 \(0\) 和 \(1\) 之間的等間距單調序列（即 均勻參數化）。

輸出包含兩個物件：一個 3 元組 \(\left(t,c,k\right)\)，其中包含樣條表示和參數變數 \(u.\)

係數是 \(N\) 個陣列的清單，其中每個陣列對應於輸入資料的一個維度。請注意，節點 t 對應於曲線 u 的參數化。

關鍵字引數 s 用於指定在樣條擬合期間要執行的平滑量。\(s\) 的預設值為 \(s=m-\sqrt{2m}\)，其中 \(m\) 是要擬合的資料點的數量。因此，如果不需要平滑，則應將值 \(\mathbf{s}=0\) 傳遞給例程。

一旦確定了資料的樣條表示，就可以使用 splev 函數或透過將 tck 元組包裝到 BSpline 物件中來評估它，如下所示。

我們先說明 s 參數對平滑某些合成雜訊資料的影響

>>> import numpy as np
>>> from scipy.interpolate import splrep, BSpline

產生一些雜訊資料

>>> x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/16)
>>> rng = np.random.default_rng()
>>> y =  np.sin(x) + 0.4*rng.standard_normal(size=len(x))

建構兩個具有不同 s 值的樣條。

>>> tck = splrep(x, y, s=0)
>>> tck_s = splrep(x, y, s=len(x))

並繪製它們

>>> import matplotlib.pyplot as plt
>>> xnew = np.arange(0, 9/4, 1/50) * np.pi
>>> plt.plot(xnew, np.sin(xnew), '-.', label='sin(x)')
>>> plt.plot(xnew, BSpline(*tck)(xnew), '-', label='s=0')
>>> plt.plot(xnew, BSpline(*tck_s)(xnew), '-', label=f's={len(x)}')
>>> plt.plot(x, y, 'o')
>>> plt.legend()
>>> plt.show()

我們看到 s=0 曲線遵循資料點的（隨機）波動，而 s > 0 曲線接近底層的正弦函數。另請注意，外插值根據 s 的值而發生巨大變化。

s 的預設值取決於是否提供權重，並且 splrep 和 splprep 的預設值也不同。因此，我們建議始終明確提供 s 的值。

操作樣條物件：程序性 (splXXX)

一旦資料的樣條表示法被確定，就可以使用函數來評估樣條 (splev) 及其導數 (splev、spalde) 在任何點的值，以及樣條在任意兩點之間的積分 ( splint)。此外，對於三次樣條 ( \(k=3\) ) 且具有 8 個或更多節點的情況，可以估計樣條的根 ( sproot)。以下範例將示範這些函數。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import interpolate

三次樣條

>>> x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8)
>>> y = np.sin(x)
>>> tck = interpolate.splrep(x, y, s=0)
>>> xnew = np.arange(0, 2*np.pi, np.pi/50)
>>> ynew = interpolate.splev(xnew, tck, der=0)

請注意，最後一行等同於 BSpline(*tck)(xnew)。

>>> plt.figure()
>>> plt.plot(x, y, 'x', xnew, ynew, xnew, np.sin(xnew), x, y, 'b')
>>> plt.legend(['Linear', 'Cubic Spline', 'True'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('Cubic-spline interpolation')
>>> plt.show()

樣條的導數

>>> yder = interpolate.splev(xnew, tck, der=1)   # or BSpline(*tck)(xnew, 1)
>>> plt.figure()
>>> plt.plot(xnew, yder, xnew, np.cos(xnew),'--')
>>> plt.legend(['Cubic Spline', 'True'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('Derivative estimation from spline')
>>> plt.show()

樣條的所有導數

>>> yders = interpolate.spalde(xnew, tck)
>>> plt.figure()
>>> for i in range(len(yders[0])):
...    plt.plot(xnew, [d[i] for d in yders], '--', label=f"{i} derivative")
>>> plt.legend()
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('All derivatives of a B-spline')
>>> plt.show()

樣條的積分

>>> def integ(x, tck, constant=-1):
...     x = np.atleast_1d(x)
...     out = np.zeros(x.shape, dtype=x.dtype)
...     for n in range(len(out)):
...         out[n] = interpolate.splint(0, x[n], tck)
...     out += constant
...     return out

>>> yint = integ(xnew, tck)
>>> plt.figure()
>>> plt.plot(xnew, yint, xnew, -np.cos(xnew), '--')
>>> plt.legend(['Cubic Spline', 'True'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('Integral estimation from spline')
>>> plt.show()

樣條的根

>>> interpolate.sproot(tck)
array([3.1416])  # may vary

請注意，sproot 可能無法在近似區間的邊緣找到明顯的解，\(x = 0\)。如果我們在稍微更大的區間上定義樣條，我們將恢復兩個根 \(x = 0\) 和 \(x = \pi\)

>>> x = np.linspace(-np.pi/4, np.pi + np.pi/4, 51)
>>> y = np.sin(x)
>>> tck = interpolate.splrep(x, y, s=0)
>>> interpolate.sproot(tck)
array([0., 3.1416])

參數樣條

>>> t = np.arange(0, 1.1, .1)
>>> x = np.sin(2*np.pi*t)
>>> y = np.cos(2*np.pi*t)
>>> tck, u = interpolate.splprep([x, y], s=0)
>>> unew = np.arange(0, 1.01, 0.01)
>>> out = interpolate.splev(unew, tck)
>>> plt.figure()
>>> plt.plot(x, y, 'x', out[0], out[1], np.sin(2*np.pi*unew), np.cos(2*np.pi*unew), x, y, 'b')
>>> plt.legend(['Linear', 'Cubic Spline', 'True'])
>>> plt.axis([-1.05, 1.05, -1.05, 1.05])
>>> plt.title('Spline of parametrically-defined curve')
>>> plt.show()

請注意，在最後一個範例中，splprep 將樣條係數作為陣列列表回傳，其中每個陣列對應於輸入資料的一個維度。因此，為了將其輸出包裝到 BSpline 中，我們需要轉置係數 (或使用 BSpline(..., axis=1))

>>> tt, cc, k = tck
>>> cc = np.array(cc)
>>> bspl = BSpline(tt, cc.T, k)    # note the transpose
>>> xy = bspl(u)
>>> xx, yy = xy.T   # transpose to unpack into a pair of arrays
>>> np.allclose(x, xx)
True
>>> np.allclose(y, yy)
True

物件導向 (UnivariateSpline)

上述描述的樣條擬合功能也可以透過物件導向介面使用。一維樣條是 UnivariateSpline 類別的物件，並使用曲線的 \(x\) 和 \(y\) 分量作為建構子的引數來建立。該類別定義了 __call__，允許使用 x 軸值呼叫物件，在這些 x 軸值上應評估樣條，並傳回內插的 y 值。以下範例針對子類別 InterpolatedUnivariateSpline 進行了展示。integral、derivatives 和 roots 方法也可在 UnivariateSpline 物件上使用，允許計算樣條的定積分、導數和根。

UnivariateSpline 類別也可以透過提供非零的平滑參數 s 來平滑資料，其含義與上述 splrep 函數的 s 關鍵字相同。這會產生一個節點數量少於資料點數量的樣條，因此不再是嚴格意義上的內插樣條，而是平滑樣條。如果不需要這樣，可以使用 InterpolatedUnivariateSpline 類別。它是 UnivariateSpline 的子類別，它總是通過所有點（相當於強制平滑參數為 0）。以下範例將示範此類別。

LSQUnivariateSpline 類別是 UnivariateSpline 的另一個子類別。它允許使用者使用參數 t 顯式指定內部節點的數量和位置。這允許建立具有非線性間距的自訂樣條，以便在某些域中進行內插，在其他域中進行平滑，或更改樣條的特性。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import interpolate

內插 UnvariateSpline

>>> x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8)
>>> y = np.sin(x)
>>> s = interpolate.InterpolatedUnivariateSpline(x, y)
>>> xnew = np.arange(0, 2*np.pi, np.pi/50)
>>> ynew = s(xnew)

>>> plt.figure()
>>> plt.plot(x, y, 'x', xnew, ynew, xnew, np.sin(xnew), x, y, 'b')
>>> plt.legend(['Linear', 'InterpolatedUnivariateSpline', 'True'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('InterpolatedUnivariateSpline')
>>> plt.show()

具有非均勻節點的 LSQUnivarateSpline

>>> t = [np.pi/2-.1, np.pi/2+.1, 3*np.pi/2-.1, 3*np.pi/2+.1]
>>> s = interpolate.LSQUnivariateSpline(x, y, t, k=2)
>>> ynew = s(xnew)

>>> plt.figure()
>>> plt.plot(x, y, 'x', xnew, ynew, xnew, np.sin(xnew), x, y, 'b')
>>> plt.legend(['Linear', 'LSQUnivariateSpline', 'True'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('Spline with Specified Interior Knots')
>>> plt.show()

二維平滑樣條

除了平滑一維樣條外，FITPACK 函式庫還提供了將二維曲面擬合到二維資料的方法。這些曲面可以被視為兩個引數的函數，\(z = g(x, y)\)，建構為一維樣條的張量積。

假設資料保存在三個陣列 x、y 和 z 中，則可以透過兩種方式解釋這些資料陣列。首先是散佈內插問題—資料被假定為成對的，即值對 x[i] 和 y[i] 代表點 i 的座標，該點對應於 z[i]。

建構曲面 \(g(x, y)\) 以滿足

\[\sum_i \left[ w_i (g(x_i, y_i) - z_i)\right]^2 \leqslant s\]

其中 \(w_i\) 是非負權重，而 s 是輸入參數，稱為平滑因子，它控制結果函數 g(x, y) 的平滑度和資料近似的品質（即 \(g(x_i, y_i)\) 和 \(z_i\) 之間的差異）之間的相互作用。\(s = 0\) 的極限在形式上對應於內插，其中曲面通過輸入資料，\(g(x_i, y_i) = z_i\)。但是，請參閱下面的註解。

第二種情況—矩形網格內插問題—是資料點被假定位於由 x 和 y 陣列的元素的所有對組定義的矩形網格上。對於此問題，z 陣列被假定為二維的，並且 z[i, j] 對應於 (x[i], y[j])。建構雙變數樣條函數 \(g(x, y)\) 以滿足

\[\sum_i \sum_j \left[ (g(x_i, y_j) - z_{i,j})\right]^2 \leqslant s\]

其中 s 是平滑因子。在這裡，\(s=0\) 的極限在形式上也對應於內插，\(g(x_i, y_j) = z_{i, j}\)。

注意

在內部，平滑曲面 \(g(x, y)\) 是透過將樣條節點放置到由資料陣列定義的邊界框中來建構的。節點是透過 FITPACK 演算法自動放置的，直到達到所需的平滑度。

節點可以放置在遠離資料點的位置。

雖然 \(s=0\) 在形式上對應於雙變數樣條內插，但 FITPACK 演算法並非旨在用於內插，並且可能導致意外的結果。

對於散佈資料內插，請優先使用 griddata；對於規則網格上的資料，請優先使用 RegularGridInterpolator。

注意

如果輸入資料 x 和 y 使得輸入維度具有不可通約的單位並且相差許多數量級，則內插器 \(g(x, y)\) 可能具有數值假影。考慮在內插之前重新縮放資料。

現在我們依次考慮兩個樣條擬合問題。

散佈資料的雙變數樣條擬合

底層 FITPACK 函式庫有兩個介面，一個程序式介面和一個物件導向介面。

程序式介面 (`bisplrep`)

對於二維曲面的（平滑）樣條擬合，可以使用函數 bisplrep。此函數將一維陣列 x、y 和 z 作為必要輸入，它們表示曲面 \(z=f(x, y).\) 上的點。可以在 x 和 y 方向上透過可選參數 kx 和 ky 指定樣條階數。預設值是雙三次樣條，kx=ky=3。

bisplrep 的輸出是一個列表 [tx ,ty, c, kx, ky]，其條目分別表示節點位置的分量、樣條的係數以及每個座標中樣條的階數。將此列表保存在單個物件 tck 中很方便，這樣可以輕鬆地將其傳遞給函數 bisplev。關鍵字 s 可用於在確定適當的樣條時更改對資料執行的平滑量。 \(s\) 的建議值取決於權重 \(w_i\)。如果這些權重取為 \(1/d_i\)，其中 \(d_i\) 是 \(z_i\) 標準差的估計值，則應該在範圍 \(m- \sqrt{2m}, m + \sqrt{2m}\) 內找到 \(s\) 的良好值，其中 \(m\) 是 x、y 和 z 向量中的資料點數。

預設值為 \(s=m-\sqrt{2m}\)。因此，如果不需要平滑，則應將 ``s=0`` 傳遞給 `bisplrep`。（但請參閱上面的註解）。

為了評估二維樣條及其偏導數（直到樣條的階數），需要函數 bisplev。此函數將兩個一維陣列作為前兩個引數，它們的叉積指定了評估樣條的域。第三個引數是從 bisplrep 回傳的 tck 列表。如果需要，第四個和第五個引數分別提供 \(x\) 和 \(y\) 方向上的偏導數階數。

注意

重要的是要注意，二維內插不應用於尋找影像的樣條表示法。所使用的演算法不適用於大量輸入點。scipy.signal 和 scipy.ndimage 包含更適合尋找影像樣條表示法的演算法。

二維內插命令旨在用於內插二維函數，如下面的範例所示。此範例使用 NumPy 中的 mgrid 命令，該命令對於在多個維度中定義「網格」非常有用。（如果不需要完整網格，另請參閱 ogrid 命令）。輸出引數的數量和每個引數的維度數由傳遞到 mgrid 中的索引物件數量決定。

>>> import numpy as np
>>> from scipy import interpolate
>>> import matplotlib.pyplot as plt

在稀疏的 20x20 網格上定義函數

>>> x_edges, y_edges = np.mgrid[-1:1:21j, -1:1:21j]
>>> x = x_edges[:-1, :-1] + np.diff(x_edges[:2, 0])[0] / 2.
>>> y = y_edges[:-1, :-1] + np.diff(y_edges[0, :2])[0] / 2.
>>> z = (x+y) * np.exp(-6.0*(x*x+y*y))

>>> plt.figure()
>>> lims = dict(cmap='RdBu_r', vmin=-0.25, vmax=0.25)
>>> plt.pcolormesh(x_edges, y_edges, z, shading='flat', **lims)
>>> plt.colorbar()
>>> plt.title("Sparsely sampled function.")
>>> plt.show()

在新的 70x70 網格上內插函數

>>> xnew_edges, ynew_edges = np.mgrid[-1:1:71j, -1:1:71j]
>>> xnew = xnew_edges[:-1, :-1] + np.diff(xnew_edges[:2, 0])[0] / 2.
>>> ynew = ynew_edges[:-1, :-1] + np.diff(ynew_edges[0, :2])[0] / 2.
>>> tck = interpolate.bisplrep(x, y, z, s=0)
>>> znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)

>>> plt.figure()
>>> plt.pcolormesh(xnew_edges, ynew_edges, znew, shading='flat', **lims)
>>> plt.colorbar()
>>> plt.title("Interpolated function.")
>>> plt.show()

物件導向介面 (`SmoothBivariateSpline`)

用於散佈資料的雙變數樣條平滑的物件導向介面，SmoothBivariateSpline 類別，實作了 bisplrep / bisplev 對組功能的一個子集，並且具有不同的預設值。

它採用權重陣列的元素等於 1，\(w_i = 1\)，並根據平滑因子 s 的輸入值自動建構節點向量—預設值為 \(m\)，即資料點的數量。

x 和 y 方向上的樣條階數由可選參數 kx 和 ky 控制，預設值為 kx=ky=3。

我們使用以下範例來說明平滑因子的效果

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import SmoothBivariateSpline

import warnings
warnings.simplefilter('ignore')

train_x, train_y = np.meshgrid(np.arange(-5, 5, 0.5), np.arange(-5, 5, 0.5))
train_x = train_x.flatten()
train_y = train_y.flatten()

def z_func(x, y):
    return np.cos(x) + np.sin(y) ** 2 + 0.05 * x + 0.1 * y

train_z = z_func(train_x, train_y)
interp_func = SmoothBivariateSpline(train_x, train_y, train_z, s=0.0)
smth_func = SmoothBivariateSpline(train_x, train_y, train_z)

test_x = np.arange(-9, 9, 0.01)
test_y = np.arange(-9, 9, 0.01)
grid_x, grid_y = np.meshgrid(test_x, test_y)

interp_result = interp_func(test_x, test_y).T
smth_result = smth_func(test_x, test_y).T
perfect_result = z_func(grid_x, grid_y)

fig, axes = plt.subplots(1, 3, figsize=(16, 8))
extent = [test_x[0], test_x[-1], test_y[0], test_y[-1]]
opts = dict(aspect='equal', cmap='nipy_spectral', extent=extent, vmin=-1.5, vmax=2.5)

im = axes[0].imshow(perfect_result, **opts)
fig.colorbar(im, ax=axes[0], orientation='horizontal')
axes[0].plot(train_x, train_y, 'w.')
axes[0].set_title('Perfect result, sampled function', fontsize=21)

im = axes[1].imshow(smth_result, **opts)
axes[1].plot(train_x, train_y, 'w.')
fig.colorbar(im, ax=axes[1], orientation='horizontal')
axes[1].set_title('s=default', fontsize=21)

im = axes[2].imshow(interp_result, **opts)
fig.colorbar(im, ax=axes[2], orientation='horizontal')
axes[2].plot(train_x, train_y, 'w.')
axes[2].set_title('s=0', fontsize=21)

plt.tight_layout()
plt.show()

在這裡，我們採用一個已知的函數（顯示在最左側的面板中），在點網格上對其進行取樣（以白色點顯示），並使用預設平滑（中間面板）和強制內插（最右側面板）建構樣條擬合。

幾個特徵清晰可見。首先，s 的預設值為此資料提供了過多的平滑；強制內插條件 s = 0 允許將底層函數恢復到合理的準確度。其次，在內插範圍之外（即，白色點覆蓋的區域），結果是使用最近鄰常數進行外插。最後，我們不得不消除警告（這是一個糟糕的形式，是的！）。

此處的警告在 s=0 的情況下發出，並表示當我們強制內插條件時，FITPACK 遇到的內部困難。如果您在程式碼中看到此警告，請考慮切換到 bisplrep 並增加其 nxest、nyest 參數（有關更多詳細資訊，請參閱 bisplrep 文件字串）。

網格上資料的雙變數樣條擬合

對於網格化的二維資料，可以使用 RectBivariateSpline 類別完成擬合平滑張量積樣條。它具有與 SmoothBivariateSpline 類似的介面，主要區別在於一維輸入陣列 x 和 y 被理解為定義二維網格（作為其外積），並且 z 陣列是二維的，其形狀為 len(x) 乘以 len(y)。

x 和 y 方向上的樣條階數由可選參數 kx 和 ky 控制，預設值為 kx=ky=3，即雙三次樣條。

平滑因子的預設值為 s=0。儘管如此，我們仍然建議始終顯式指定 s。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import RectBivariateSpline

x = np.arange(-5.01, 5.01, 0.25)        # the grid is an outer product
y = np.arange(-5.01, 7.51, 0.25)        # of x and y arrays

xx, yy = np.meshgrid(x, y, indexing='ij')
z = np.sin(xx**2 + 2.*yy**2)            # z array needs to be 2-D

func = RectBivariateSpline(x, y, z, s=0)

xnew = np.arange(-5.01, 5.01, 1e-2)
ynew = np.arange(-5.01, 7.51, 1e-2)
znew = func(xnew, ynew)

plt.imshow(znew)
plt.colorbar()
plt.show()

球座標中資料的雙變數樣條擬合

如果您的資料以球座標給出，\(r = r(\theta, \phi)\)，SmoothSphereBivariateSpline 和 RectSphereBivariateSpline 分別提供了 SmoothBivariateSpline 和 RectBivariateSpline 的方便類比。

這些類別確保樣條擬合對於 \(\theta \in [0, \pi]\) 和 \(\phi \in [0, 2\pi]\) 的週期性，並提供對極點連續性的一些控制。有關詳細資訊，請參閱這些類別的文件字串。