Fitting a monotonically increasing spline function

Question

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

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

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

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?

Answer 1

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

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)