python 相关性分析切点寻找,Python自然平滑样条线

I am trying to find a python package that would give an option to fit natural smoothing splines with user selectable smoothing factor. Is there an implementation for that? If not, how would you use what is available to implement it yourself?

By natural spline I mean that there should be a condition that the second derivative of the fitted function at the endpoints is zero (linear).

By smoothing spline I mean that the spline should not be 'interpolating' (passing through all the datapoints). I would like to decide the correct smoothing factor lambda (see the Wikipedia page for smoothing splines) myself.

What I have found

scipy.interpolate.CubicSpline [link]: Does natural (cubic) spline fitting. Does interpolation, and there is no way to smooth the data.

scipy.interpolate.UnivariateSpline [link]: Does spline fitting with user selectable smoothing factor. However, there is no option to make the splines natural.

解决方案

After hours of investigation, I did not find any pip installable packages which could fit a natural cubic spline with user-controllable smoothness. However, after deciding to write one myself, while reading about the topic I stumbled upon a blog post by github user madrury. He has written python code capable of producing natural cubic spline models.

The model code is available here (NaturalCubicSpline) with a BSD-licence. He has also written some examples in an IPython notebook.

But since this is the Internet and links tend to die, I will copy the relevant parts of the source code here + a helper function (get_natural_cubic_spline_model) written by me, and show an example of how to use it. The smoothness of the fit can be controlled by using different number of knots. The position of the knots can be also specified by the user.

Example

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()

The source code for 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