Python自然平滑样条线

Question

我试图找到一个python软件包，该软件包将提供一个选项，以使自然平滑样条线与用户可选的平滑因子相匹配。有没有实现的方法？如果没有，您将如何使用可用的工具自己实施？

• 所谓自然样条曲线，是指应该满足以下条件：拟合函数在端点处的二阶导数为零（线性）。

• 通过平滑样条曲线，我的意思是样条曲线不应被“插值”（通过所有数据点）。我想自己决定正确的平滑系数lambda（请参见Wikipedia页面以平滑样条线）。

我发现了什么

• scipy.interpolate.CubicSpline [ 链接 ]：执行自然（三次）样条曲线拟合。进行插值，并且无法平滑数据。

• scipy.interpolate.UnivariateSpline [ 链接 ]：样条曲线是否与用户可选的平滑因子拟合。但是，无法选择使花键自然化。

Answer 1

经过数小时的调查，我没有找到任何可安装在用户喜欢的平滑度的自然三次样条上的可点子安装的软件包。但是，在决定写一个自己的文章后，在阅读有关该主题的文章时，我偶然发现了github用户madrury的博客文章。他编写了能够生成自然三次样条模型的python代码。

此处具有BSD许可证的模型代码（NaturalCubicSpline）。他还在IPython笔记本中编写了一些示例。

但是由于这是Internet，并且链接容易消失，因此我将在此处复制源代码的相关部分+ 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

Answer 2

python 包 patsy 具有生成样条基的函数，包括自然三次样条基。文档中描述。然后可以使用任何库来拟合模型，例如 scikit-learn 或 statsmodels。

• 该df参数cr()可用于控制“平滑度”
• 请注意，太低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()

Answer 3

您可以将这种自然三次平滑样条的numpy / scipy实现用于单变量/多变量数据平滑。平滑参数应在[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()

Answer 4

对于我的一个项目，我需要为时间序列建模创建间隔，并且为了使程序更高效，我创建了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 可以以矢量化的方式对多个时间序列进行平滑处理