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 准则来找到最佳正则化参数,否则我们使用给定的参数解决正则化加权线性回归问题。该参数以以下方式控制权衡:参数越大,函数越平滑。- 参数:
- xarray_like,形状 (n,)
横坐标。n 必须至少为 5。
- yarray_like,形状 (n,)
纵坐标。n 必须至少为 5。
- warray_like,形状 (n,),可选
权重向量。默认值为
np.ones_like(x)
。- lamfloat,(\(\lambda \geq 0\)), 可选
正则化参数。如果
lam
为 None,则从 GCV 准则中找到它。默认值为 None。
- 返回:
- func一个 BSpline 对象。
一个可调用的对象,表示 B 样条基中的样条,作为使用 GCV 准则 [1] 平滑样条问题的解决方案,如果
lam
为 None,否则使用给定的参数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
中,该和由权重的平方构建。如果初始问题不适定(例如,\(X^T W X\) 的乘积,其中 \(X\) 是设计矩阵,不是正定矩阵),则会引发 ValueError。
参考文献
[1]G. Wahba,“估计平滑参数”,载于《观测数据的样条模型》,费城,宾夕法尼亚州:工业与应用数学学会,1990 年,第 45-65 页。DOI:10.1137/1.9781611970128
[2]H. J. Woltring,用于广义交叉验证样条平滑和微分的 Fortran 包,《工程软件进展》,第 8 卷,第 2 期,第 104-113 页,1986 年。DOI:10.1016/0141-1195(86)90098-7
[3]T. Hastie、J. Friedman 和 R. Tisbshirani,“平滑样条”,载于《统计学习要素:数据挖掘、推断和预测》,纽约:施普林格,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()