Python 自然平滑样条曲线
Posted
技术标签:
【中文标题】Python 自然平滑样条曲线【英文标题】：Python natural smoothing splines 【发布时间】：2018-12-21 14:20:18 【问题描述】：我正在尝试找到一个 python 包，它可以提供一个选项来适应 自然平滑样条 与用户可选择的平滑因子。有没有实现呢？如果没有，您将如何使用可用的东西来自己实现它？
通过自然样条，我的意思是应该有一个条件，即在端点处拟合函数的二阶导数为零（线性）。
通过平滑样条曲线我的意思是样条曲线不应该是“插值”（通过所有数据点）。我想自己决定正确的平滑因子 lambda（参见 Wikipedia page 平滑样条曲线）。
我发现了什么
scipy.interpolate.CubicSpline [link]：自然（三次）样条拟合。做插值，没办法平滑数据。
scipy.interpolate.UnivariateSpline [link]：样条曲线是否与用户可选择的平滑因子拟合。但是，没有使样条线自然的选项。
【问题讨论】：
【参考方案1】：python 包 patsy 具有生成样条基的功能，包括自然三次样条基。在documentation 中描述。 然后可以使用任何库来拟合模型，例如scikit-learn 或 statsmodels。
cr() 的df 参数可用于控制“平滑度”
请注意，太低的df 可能会导致拟合不足（见下文）。
一个使用 scikit-learn 的简单示例。
import numpy as np
from sklearn.linear_model import LinearRegression
from patsy import cr
import matplotlib.pyplot as plt
n_obs = 600
np.random.seed(0)
x = np.linspace(-3, 3, n_obs)
y = 1 / (x ** 2 + 1) * np.cos(np.pi * x) + np.random.normal(0, 0.2, size=n_obs)
def plot_smoothed(df=5):
# Generate spline basis with different degrees of freedom
x_basis = cr(x, df=df, constraints="center")
# Fit model to the data
model = LinearRegression().fit(x_basis, y)
# Get estimates
y_hat = model.predict(x_basis)
plt.plot(x, y_hat, label=f"df=df")
plt.scatter(x, y, s=4, color="tab:blue")
for df in (5, 7, 10, 25):
plot_smoothed(df)
plt.legend()
plt.title(f"Natural cubic spline with varying degrees of freedom")
plt.show()
【讨论】：
【参考方案2】：对于我的一个项目，我需要为时间序列建模创建间隔，并且为了使过程更高效，我创建了tsmoothie：一个用于以矢量化方式进行时间序列平滑和异常值检测的 python 库。
它提供了不同的平滑算法以及计算间隔的可能性。
对于自然立方类型的SplineSmoother：
import numpy as np
import matplotlib.pyplot as plt
from tsmoothie.smoother import *
def func(x):
return 1/(1+25*x**2)
# make example data
x = np.linspace(-1,1,300)
y = func(x) + np.random.normal(0, 0.2, len(x))
# operate smoothing
smoother = SplineSmoother(n_knots=10, spline_type='natural_cubic_spline')
smoother.smooth(y)
# generate intervals
low, up = smoother.get_intervals('prediction_interval', confidence=0.05)
# plot the first smoothed timeseries with intervals
plt.figure(figsize=(11,6))
plt.plot(smoother.smooth_data[0], linewidth=3, color='blue')
plt.plot(smoother.data[0], '.k')
plt.fill_between(range(len(smoother.data[0])), low[0], up[0], alpha=0.3)
我还指出tsmoothie可以以矢量化的方式对多个时间序列进行平滑
【讨论】：
【参考方案3】：编程语言 R 提供了一个非常好的自然三次平滑样条实现。您可以通过rpy2 在 Python 中使用 R 函数：
import rpy2.robjects as robjects
r_y = robjects.FloatVector(y_train)
r_x = robjects.FloatVector(x_train)
r_smooth_spline = robjects.r['smooth.spline'] #extract R function# run smoothing function
spline1 = r_smooth_spline(x=r_x, y=r_y, spar=0.7)
ySpline=np.array(robjects.r['predict'](spline1,robjects.FloatVector(x_smooth)).rx2('y'))
plt.plot(x_smooth,ySpline)
如果要直接设置lambda：spline1 = r_smooth_spline(x=r_x, y=r_y, lambda=42)是不行的，因为lambda在Python中已经有别的意思了，但是有一个解决办法：How to use the lambda argument of smooth.spline in RPy WITHOUT Python interprating it as lambda。
要运行代码，您首先需要定义数据x_train 和y_train，如果您想以全高清分辨率在-3 和5 之间绘制它，您可以定义x_smooth=np.array(np.linspace(-3,5,1920)).
。
【讨论】：
【参考方案4】：您可以使用自然三次平滑样条的this numpy/scipy implementation 进行单变量/多变量数据平滑。平滑参数应在 [0.0, 1.0] 范围内。如果我们使用等于 1.0 的平滑参数，我们会得到自然三次样条插值，而无需数据平滑。该实现还支持单变量数据的矢量化。
单变量示例：
import numpy as np
import matplotlib.pyplot as plt
import csaps
np.random.seed(1234)
x = np.linspace(-5., 5., 25)
y = np.exp(-(x/2.5)**2) + (np.random.rand(25) - 0.2) * 0.3
sp = csaps.UnivariateCubicSmoothingSpline(x, y, smooth=0.85)
xs = np.linspace(x[0], x[-1], 150)
ys = sp(xs)
plt.plot(x, y, 'o', xs, ys, '-')
plt.show()
双变量示例：
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import csaps
xdata = [np.linspace(-3, 3, 61), np.linspace(-3.5, 3.5, 51)]
i, j = np.meshgrid(*xdata, indexing='ij')
ydata = (3 * (1 - j)**2. * np.exp(-(j**2) - (i + 1)**2)
- 10 * (j / 5 - j**3 - i**5) * np.exp(-j**2 - i**2)
- 1 / 3 * np.exp(-(j + 1)**2 - i**2))
np.random.seed(12345)
noisy = ydata + (np.random.randn(*ydata.shape) * 0.75)
sp = csaps.MultivariateCubicSmoothingSpline(xdata, noisy, smooth=0.988)
ysmth = sp(xdata)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(j, i, noisy, linewidths=0.5, color='r')
ax.scatter(j, i, noisy, s=5, c='r')
ax.plot_surface(j, i, ysmth, linewidth=0, alpha=1.0)
plt.show()
【讨论】：
csaps 看起来很棒！需要注意的一件事是，至少在 v. 0.11.0 中，csaps.csaps() 需要唯一的 x 值。我最终有一个预处理步骤来计算非唯一 x 值的平均值。
@np8，x 值必须是唯一且严格单调的。如果您有非唯一或/和非单调的 x 值，则需要使用带有参数化的多元平滑 xs = u(x), ys = v(y)。
我收到此错误：AttributeError: module 'csaps' has no attribute 'MultivariateCubicSmoothingSpline' 我找不到原因！
@sam，请阅读当前文档。 csaps.readthedocs.io/en/latest
@iroln 感谢您的评论。我仍然无法解决问题。在文档中没有多变量案例的代码或示例（只有一个数字！）。将 import csaps 更改为 from csaps import csaps 也没有帮助。【参考方案5】：
经过数小时的调查，我没有找到任何可以安装自然三次样条的 pip 可安装包，并且具有用户可控的平滑度。然而，在决定自己写一个之后，在阅读该主题时，我偶然发现了 github 用户 madrury 的 blog post。他编写了能够生成自然三次样条模型的 Python 代码。
型号代码在here (NaturalCubicSpline) 和a BSD-licence 可用。他还在IPython notebook 中写了一些示例。
但是由于这是互联网，链接容易死，我将这里的源代码的相关部分+我编写的辅助函数（get_natural_cubic_spline_model），并显示如何使用它的示例。可以通过使用不同数量的结来控制拟合的平滑度。结的位置也可以由用户指定。
示例
from matplotlib import pyplot as plt
import numpy as np
def func(x):
return 1/(1+25*x**2)
# make example data
x = np.linspace(-1,1,300)
y = func(x) + np.random.normal(0, 0.2, len(x))
# The number of knots can be used to control the amount of smoothness
model_6 = get_natural_cubic_spline_model(x, y, minval=min(x), maxval=max(x), n_knots=6)
model_15 = get_natural_cubic_spline_model(x, y, minval=min(x), maxval=max(x), n_knots=15)
y_est_6 = model_6.predict(x)
y_est_15 = model_15.predict(x)
plt.plot(x, y, ls='', marker='.', label='originals')
plt.plot(x, y_est_6, marker='.', label='n_knots = 6')
plt.plot(x, y_est_15, marker='.', label='n_knots = 15')
plt.legend(); plt.show()
get_natural_cubic_spline_model的源代码
import numpy as np
import pandas as pd
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
def get_natural_cubic_spline_model(x, y, minval=None, maxval=None, n_knots=None, knots=None):
"""
Get a natural cubic spline model for the data.
For the knots, give (a) `knots` (as an array) or (b) minval, maxval and n_knots.
If the knots are not directly specified, the resulting knots are equally
space within the *interior* of (max, min). That is, the endpoints are
*not* included as knots.
Parameters
----------
x: np.array of float
The input data
y: np.array of float
The outpur data
minval: float
Minimum of interval containing the knots.
maxval: float
Maximum of the interval containing the knots.
n_knots: positive integer
The number of knots to create.
knots: array or list of floats
The knots.
Returns
--------
model: a model object
The returned model will have following method:
- predict(x):
x is a numpy array. This will return the predicted y-values.
"""
if knots:
spline = NaturalCubicSpline(knots=knots)
else:
spline = NaturalCubicSpline(max=maxval, min=minval, n_knots=n_knots)
p = Pipeline([
('nat_cubic', spline),
('regression', LinearRegression(fit_intercept=True))
])
p.fit(x, y)
return p
class AbstractSpline(BaseEstimator, TransformerMixin):
"""Base class for all spline basis expansions."""
def __init__(self, max=None, min=None, n_knots=None, n_params=None, knots=None):
if knots is None:
if not n_knots:
n_knots = self._compute_n_knots(n_params)
knots = np.linspace(min, max, num=(n_knots + 2))[1:-1]
max, min = np.max(knots), np.min(knots)
self.knots = np.asarray(knots)
@property
def n_knots(self):
return len(self.knots)
def fit(self, *args, **kwargs):
return self
class NaturalCubicSpline(AbstractSpline):
"""Apply a natural cubic basis expansion to an array.
The features created with this basis expansion can be used to fit a
piecewise cubic function under the constraint that the fitted curve is
linear *outside* the range of the knots.. The fitted curve is continuously
differentiable to the second order at all of the knots.
This transformer can be created in two ways:
- By specifying the maximum, minimum, and number of knots.
- By specifying the cutpoints directly.
If the knots are not directly specified, the resulting knots are equally
space within the *interior* of (max, min). That is, the endpoints are
*not* included as knots.
Parameters
----------
min: float
Minimum of interval containing the knots.
max: float
Maximum of the interval containing the knots.
n_knots: positive integer
The number of knots to create.
knots: array or list of floats
The knots.
"""
def _compute_n_knots(self, n_params):
return n_params
@property
def n_params(self):
return self.n_knots - 1
def transform(self, X, **transform_params):
X_spl = self._transform_array(X)
if isinstance(X, pd.Series):
col_names = self._make_names(X)
X_spl = pd.DataFrame(X_spl, columns=col_names, index=X.index)
return X_spl
def _make_names(self, X):
first_name = "_spline_linear".format(X.name)
rest_names = ["_spline_".format(X.name, idx)
for idx in range(self.n_knots - 2)]
return [first_name] + rest_names
def _transform_array(self, X, **transform_params):
X = X.squeeze()
try:
X_spl = np.zeros((X.shape[0], self.n_knots - 1))
except IndexError: # For arrays with only one element
X_spl = np.zeros((1, self.n_knots - 1))
X_spl[:, 0] = X.squeeze()
def d(knot_idx, x):
def ppart(t): return np.maximum(0, t)
def cube(t): return t*t*t
numerator = (cube(ppart(x - self.knots[knot_idx]))
- cube(ppart(x - self.knots[self.n_knots - 1])))
denominator = self.knots[self.n_knots - 1] - self.knots[knot_idx]
return numerator / denominator
for i in range(0, self.n_knots - 2):
X_spl[:, i+1] = (d(i, X) - d(self.n_knots - 2, X)).squeeze()
return X_spl
【讨论】：
@np8 这是仅考虑训练数据来拟合模型的准确方法吗？ model_6 = get_natural_cubic_spline_model(X_train, y_train, minval=min(X_train), maxval=max(X_train), n_knots=6)，其中 X_train 和 y_train 是 X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, shuffle=False)以上是关于Python 自然平滑样条曲线的主要内容，如果未能解决你的问题，请参考以下文章
是否可以在 gnuplot 中的曲线的两端使样条曲线连续且平滑
R语言广义加性模型GAMs：可视化每个变量的样条函数样条函数与变量与目标变量之间的平滑曲线比较并进行多变量的归一化比较测试广义线性加性模型GAMs在测试集上的表现(防止过拟合)