Python SciPy interpolate.make_smoothing_spline用法及代碼示例

本文簡要介紹 python 語言中 scipy.interpolate.make_smoothing_spline 的用法。

用法:

scipy.interpolate.make_smoothing_spline(x, y, w=None, lam=None)#

使用 lam 計算平滑三次樣條函數(的係數)，以控製曲線的平滑度與其與數據的接近程度之間的權衡。如果 lam 為 None，則使用 GCV 標準 [1] 來查找它。

找到平滑樣條作為正則化加權線性回歸問題的解：

\[\sum\limits_{i=1}^n w_i\lvert y_i - f(x_i) \rvert^2 + \lambda\int\limits_{x_1}^{x_n} (f^{(2)}(u))^2 d u\]

其中\(f\) 是樣條函數，\(w\) 是權重向量，\(\lambda\) 是正則化參數。

如果lam為None，我們使用GCV標準來尋找最佳正則化參數，否則我們用給定參數求解正則化加權線性回歸問題。參數通過以下方式控製權衡：參數越大，函數越平滑。

參數 ：：

x： 數組, 形狀 (n,)

橫坐標。 n 必須大於 5。

y： 數組, 形狀 (n,)

縱坐標。 n 必須大於 5。

w： 數組，形狀(n，)，可選

權重向量。默認為 np.ones_like(x) 。

lam： 浮點數，(\(\lambda \geq 0\))，可選

正則化參數。如果lam為None，則從GCV標準中找到。默認為“無”。

返回 ：：

func： BSpline 對象。

表示基於 B-spline 的樣條線的可調用函數，作為在 lam 為 None 的情況下使用 GCV 標準 [1] 平滑樣條線問題的解決方案，否則使用給定參數 lam 。

注意：

該算法是 Woltring 在 FORTRAN [2] 中引入的算法的幹淨室重新實現。由於許可證問題，原始版本無法在SciPy源代碼中使用。此處討論了重新實現的細節(僅提供俄語版本)[4]。

如果權重向量w為None，我們假設所有點的權重相等，並且權重向量是1的向量。

請注意，在加權殘差平方和中，權重不是平方：\(\sum\limits_{i=1}^n w_i\lvert y_i - f(x_i) \rvert^2\)，而在 splrep 中，總和是根據權重平方構建的。

如果初始問題是ill-posed(例如，產品\(X^T W X\)(其中\(X\) 是設計矩陣)不是正定義矩陣)，則會引發ValueError。

參考：

[1]

G. Wahba，觀測數據樣條模型中的“估計平滑參數”，賓夕法尼亞州費城：工業與應用數學學會，1990 年，第 45-65 頁。 DOI:10.1137/1.9781611970128

[2]

H. J. Woltring，用於廣義的 Fortran 包，cross-validatory 樣條平滑和微分，工程軟件進展，卷。 8、不。 2，第 104-113 頁，1986 年。DOI:10.1016/0141-1195(86)90098-7

[3]

T. Hastie、J. Friedman 和 R. Tisbshirani，“Smoothing Splines”，《統計學習的要素：數據挖掘、推理和預測》，紐約：Springer，2017 年，第 241-249 頁。 DOI:10.1007/978-0-387-84858-7

[4]

E. Zemlyanoy，“廣義交叉驗證平滑樣條”，理學學士論文，2022 年。https://www.hse.ru/ba/am/students/diplomas/620910604(俄語)

例子：

生成一些噪聲數據

>>> import numpy as np
>>> np.random.seed(1234)
>>> n = 200
>>> def func(x):
...    return x**3 + x**2 * np.sin(4 * x)
>>> x = np.sort(np.random.random_sample(n) * 4 - 2)
>>> y = func(x) + np.random.normal(scale=1.5, size=n)

製作平滑樣條函數

>>> from scipy.interpolate import make_smoothing_spline
>>> spl = make_smoothing_spline(x, y)

繪製兩者

>>> import matplotlib.pyplot as plt
>>> grid = np.linspace(x[0], x[-1], 400)
>>> plt.plot(grid, spl(grid), label='Spline')
>>> plt.plot(grid, func(grid), label='Original function')
>>> plt.scatter(x, y, marker='.')
>>> plt.legend(loc='best')
>>> plt.show()