espdev / csaps

Cubic spline approximation (smoothing)
https://csaps.readthedocs.io
MIT License
164 stars 27 forks source link

csaps.CubicSmoothingSpline does not perform natural boundary condition extrapolation correctly #52

Open shusheer opened 3 years ago

shusheer commented 3 years ago

csaps.CubicSmoothingSpline purports to perform extrapolation using the "natural" boundary condition.

From https://github.com/espdev/csaps/blob/master/docs/formulation.rst: "csaps spline is cubic only and it has natural boundary condition type."

This means the first derivative should be constant outside the range of observations provided to the smoother. This is not currently the behaviour.

Note that the scipy.interpolate.CubicSpline also exhibits the same behaviour, even when we explicitly specify bc_type='natural'

Example code to test:

import matplotlib.pyplot as plt

import numpy as np

from csaps import CubicSmoothingSpline

from scipy.interpolate import CubicSpline

xs = [1,2,3,4,5,6,7,8]
ys = [4.5, 3.6, 1.6, 0.0, -3.3, -3.1, -1.8, -1.7]

csaps_smoother = CubicSmoothingSpline(xs, ys)

scipy_smoother = CubicSpline(xs,ys, bc_type='natural')

extrapolated = np.linspace(min(xs)-10,max(xs)+10,int(21+(max(xs)-min(xs))))

fig, axs = plt.subplots(2, 2,figsize=(12,12))
axs[0, 0].plot(xs, ys,'.')
axs[0, 0].plot(extrapolated,csaps_smoother(extrapolated),'-')
axs[0, 0].set_title('csaps, deriv=0')
axs[0, 1].plot(xs, ys,'.')
axs[0, 1].plot(extrapolated,scipy_smoother(extrapolated),'-')
axs[0, 1].set_title('scipy, deriv=0')
axs[1, 0].plot(xs, csaps_smoother(xs,nu=1),'.')
axs[1, 0].plot(extrapolated,csaps_smoother(extrapolated,nu=1),'-')
axs[1, 0].set_title('csaps, deriv=1')
axs[1, 1].plot(xs, scipy_smoother(xs,nu=1),'.')
axs[1, 1].plot(extrapolated,scipy_smoother(extrapolated,nu=1),'-')
axs[1, 1].set_title('scipy, deriv=1')
plt.show()

I believe this issue is a superset of Issue #46 (periodic functions) wherein csaps.CubicSmoothingSpline._make_spline does not take account of boundary condition.

shusheer commented 3 years ago

See also https://github.com/scipy/scipy/issues/14472

ev-br commented 3 years ago

Commented in the scipy issue. Basically, CubicSpline only guarantees the value of derivatives at the boundary, and its extrapolation mode may be not very useful. However it provides all the info for a user to do the extrapolation, and it needs to be taken care of by the user.