Python SciPy的UnivariateSpline与R smooth.spline的比较

Question

我正在将一份用R语言编写的脚本移植到Python。在R中，我使用了smooth.spline方法，在Python中，我使用的是SciPy UnivariateSpline方法。尽管它们都基于三次样条方法，但它们产生的结果并不相同。是否有一种方法或替代UnivariateSpline的方法，可以使Python的样条返回与R相同的结果？

我是一名数学家。我理解样条的一般思想，但不了解Python或R中其实现的细节。

以下是用R和Python编写的代码。输入数据对两者都是相同的。

这是输入数据：

x =  0.0,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1.0
y =   -1,    1,    1,   -1,    1,    0,   .5,   .5,   .4,   .5,   -1

以下是R语言代码：

x = seq(0,1, by = .1); 
y = c(-1,1,1, -1,1,0, .5,.5,.4,  .5, -1);
spline_xy = smooth.spline(x,y)
predict(spline_xy,x)

输出结果为：

$x
 [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

$y
 [1]  0.120614583  0.170800975  0.210954680  0.238032338  0.253672155
 [6]  0.253684815  0.236432643  0.200264536  0.145403302  0.074993797
[11] -0.004853825

以下是Python代码

import numpy as np
from scipy.interpolate import UnivariateSpline
x = np.linspace(0, 1, num = 11, endpoint=True)    
y = np.array([-1,1,1, -1,1,0, .5,.5,.4,  .5, -1]) 
spline_xy = UnivariateSpline(x,y)
print('x =', x)
print('ysplined =',spline_xy(x))

输出结果为：

x = [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]

ysplined = 
[-0.26433566 -0.02587413  0.18857809 0.36585082  0.49277389  
  0.55617716 0.54289044  0.43974359  0.23356643 -0.08881119 
 -0.54055944]

我希望R中的输出结果$y和Python中的ysplined结果是一致的，但它们并不相同。希望能得到一些帮助，例如如何设置参数或者解释说明。在此提前感谢您的帮助！

- NoName

Answer 1

这些方法在我看来是不同的平滑方法。

R中的smooth.spline是一种“平滑样条”，它是一种超参数化自然样条（在每个数据点处有结点，内部为三次样条，线性外推），使用惩罚最小二乘法选择参数。您可以阅读帮助页面以了解如何计算惩罚的详细信息。

另一方面，Python的UnivariateSpline似乎根据此处的文档： https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.html是一个回归样条，通过最小二乘法拟合而成，没有惩罚项。它似乎能够自适应地选择结点数。

这些是完全不同的算法，我不会期望它们给出相等的结果。 我不知道是否有一个R软件包使用与Python相同的自适应结点选择。这个答案：https://dev59.com/q1UK5IYBdhLWcg3wwiOW#55481248声称在Python中引用了自然平滑样条实现，但我不知道它是否与R的实现匹配。

- user2554330

Answer 2

您可以使用 rpy2 在 Python 中调用 R 函数：

import numpy as np
import rpy2.robjects as robjects
x = np.linspace(0, 1, num = 11, endpoint=True)    
y = np.array([-1,1,1, -1,1,0, .5,.5,.4,  .5, -1])

r_x = robjects.FloatVector(x)
r_y = robjects.FloatVector(y)
r_smooth_spline = robjects.r['smooth.spline'] #extract R function
spline_xy = r_smooth_spline(x=r_x, y=r_y)
print('x =', x)
print('ysplined =',np.array(robjects.r['predict'](spline_xy,robjects.FloatVector(x)).rx2('y')))

输出结果如下：

x = [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
ysplined = [ 0.12061458  0.17080098  0.21095468  0.23803234  0.25367215  0.25368481
0.23643264  0.20026454  0.1454033   0.0749938  -0.00485382]

就像您所希望的那样。

_{如果您想直接设置lambda： spline_xy = r_smooth_spline(x=r_x, y=r_y, lambda=42)是无法工作的，因为lambda在Python中已经有另一种含义，但有一个解决方案：如何在RPy中使用平滑.spline的lambda参数，而不让Python将其解释为lambda。}

_{请注意，此代码与最新版本的rpy2不完全兼容Jupyter-notebooks。您可以按照NotImplementedError: Conversion 'rpy2py' not defined for objects of type '<class 'rpy2.rinterface.SexpClosure'>' only after I run the code twice中所述使用!pip install -Iv rpy2==3.4.2来修复此问题。}

- Jakob