0
$\begingroup$

I want to fit a monotonically increasing smooth spline function for a dataset

Code:

from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

x = [0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 20.85, 21.]
    y = [2.83811035, 2.81541896, 3.14311655, 3.22373554, 3.43033456, 3.50433385, 3.66794514, 3.462296, 3.59480959,
         3.56250726, 3.6209845,  3.63034523, 3.68238915, 3.69096892, 3.75560395, 3.83545191, 3.90419498]

    plt.plot(x, y, '*')
    f = interp1d(x, y, kind='cubic')
    yinp = f(x)
    plt.plot(x, yinp)
    
    plt.show()

enter image description here

    f = interp1d(x, y, kind='cubic')
    yinp = f(x)
    plt.plot(x, yinp)
    plt.show()

enter image description here

The current fit looks like the above. I would like to know how to fit a monotonically increasing spline function.

I found an example in r posted here https://stackoverflow.com/questions/25447999/how-to-make-monotonic-increasing-smooth-spline-with-smooth-spline-function. I am not sure what's the appropriate function in the scipy library.

Suggestions will be really helpful.

Crossposted: https://stackoverflow.com/questions/73605390/fitting-a-monotonously-increasing-spline-function-using-scipy

EDIT: I'm looking for something like the below (ref.)

enter image description here

EDIT 2: I could get the coeffs and knots but I am not sure how to use the coefficients and manually generate the function of the spline curve. Could someone please add a bit more detail to this?

For example, when we have 4 data points x = [0., 0.75, 1.8, 2.25] y = [2.83811035, 2.81541896, 3.14311655, 3.22373554]

I would like to print the piecewise polynomial function to understand how the spline function looks like.

EDIT 3: The solution posted below works great.

I am trying to print the spline for each segment

f0 = lambda x: p.c[0, i] * (x - p.x[i]) ** 3 + p.c[1, i] * (x - p.x[i]) ** 2 + p.c[2, i] * (x - p.x[i]) + p.c[3, i]

f0 = lambda x: [p.c[:, i] * (x - p.x[i]) ** (3 - i) for i in range(k + 1)] 
print(f0)

This prints <function fit_spline1.. at 0x0000028697B94F70>

Instead, I would like to see the cubic polynomial. Could someone please suggest how to print this out?

$\endgroup$
2
  • 1
    $\begingroup$ scipy.interpolate.pchip_interpolate is a monotonic cubic interpolator. It creates monotonic splines whenever your dataset is monotonic. But note that your dataset is not monotonic. $\endgroup$
    – Pepe
    Sep 5, 2022 at 20:12
  • $\begingroup$ @Pepe Thank you, please check my edit $\endgroup$
    – Natasha
    Sep 6, 2022 at 2:26

1 Answer 1

3
$\begingroup$

A smoothing spline might be good enough in your case. For example, scipy.interpolate.UnivariateSpline implements this. You can use it in the following way:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
import scipy 

x = np.array([0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 
              7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 
              20.85, 21.])
y = np.array([2.83811035, 2.81541896, 3.14311655, 
              3.22373554, 3.43033456, 3.50433385, 
              3.66794514, 3.462296,   3.59480959, 
              3.56250726, 3.6209845,  3.63034523, 
              3.68238915, 3.69096892, 3.75560395, 
              3.83545191, 3.90419498])
              

k=3 # poly-order
spl = UnivariateSpline(x, y, s=0.05,k=k)
xs  = np.linspace(x.min(), x.max(), 100)
plt.plot(x, y, 'ro', ms=5)
plt.plot(xs, spl(xs), 'cyan', lw=5,alpha=0.3)
plt.show()

enter image description here

By adjusting the smoothing parameter s you can adjust your fit. The downside is, that this implementation doesn't guarantee for a monotonically increasing spline function.

To get the coeffs of each spline segment in the usual power basis you can do the following (see also https://stackoverflow.com/a/61622711/13192621)

# get spline coeffs and knots
tck = (spl._data[8], spl._data[9], k)
p = scipy.interpolate.PPoly.from_spline(tck)

# plot each segment and return knots and coeffs
for i in range(k,len(spl.get_knots())+k-1):
    xs = np.linspace(p.x[i], p.x[i+1], 100)
    plt.plot(xs, np.polyval(p.c[:,i], xs - p.x[i]))
    print("knot ", p.x[i], " to ", p.x[i+1])
    print("coeffs ", p.c[:,i], "\n")

that returns

knot  0.0  to  3.75
coeffs  [ 1.49573322e-05 -2.09135260e-02  2.52383888e-01  2.75803708e+00] 

knot  3.75  to  7.5
coeffs  [ 1.85645802e-03 -2.07452560e-02  9.61634550e-02  3.41116897e+00] 

knot  7.5  to  21.0
coeffs  [6.59609159e-06 1.39896655e-04 1.88933574e-02 3.57795091e+00] 

which are the spline coeffs (highest power coeff on the left)

enter image description here

$\endgroup$
8
  • 1
    $\begingroup$ Thanks a lot, this really helps. I found an implementation here stackoverflow.com/a/28531283/8281509 which uses bézier curve. Could this guarantee a monotonically increasing spline function? $\endgroup$
    – Natasha
    Sep 6, 2022 at 5:13
  • $\begingroup$ could we get the polynomial function of the fitted curve? $\endgroup$
    – Natasha
    Sep 6, 2022 at 9:42
  • $\begingroup$ @Natasha Without any further assumptions, I don't see how a bézier curve could guarantee monotonicity. To get the spline coefficients you can return the knots by spl.get_knots() and the free coefficients of each spline segment as spl.get_coeffs(). So for each spline segment you have a cubic polynomial. $\endgroup$
    – Pepe
    Sep 6, 2022 at 15:01
  • $\begingroup$ Thank you, I could get the coeffs and knots but I am not sure how to use the coefficients and manually generate the function of the spline curve. Could you please add a bit more detail to this? $\endgroup$
    – Natasha
    Sep 7, 2022 at 13:28
  • 1
    $\begingroup$ If you don't give an argument to the lambda function f0. In this form, you just print the lambda object but don't execute it. Try f0 = lambda x, k, n, p: np.sum([p.c[i,k+n]*(x-p.x[k+n])**(k-i) for i in range(k+1)], axis=0) where x is the input range, k is the spline order, n the index of the spline segment and p the interpolate object as in the code above. Then you can get the input range for the first segment with n=0 and x = np.linspace(p.x[k+n], p.x[k+n+1], 100). To plot it you execute plt.plot(x, f0(x,k,n,p)). $\endgroup$
    – Pepe
    Sep 8, 2022 at 16:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.