# Fitting a monotonically increasing spline function

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.

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?

• 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.
– Pepe
Sep 5, 2022 at 20:12
• @Pepe Thank you, please check my edit Sep 6, 2022 at 2:26

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)

• 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? Sep 6, 2022 at 5:13
• could we get the polynomial function of the fitted curve? Sep 6, 2022 at 9:42
• @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.
– Pepe
Sep 6, 2022 at 15:01
• 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? Sep 7, 2022 at 13:28
• 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)).
– Pepe
Sep 8, 2022 at 16:04