I have some data that can be fit by eye with a linear spline with a single internal knot pretty easily. But scipy.interpolate.splrep is choking on the fit, no matter what smoothing parameter I provide. You can see a plot of my data below.
Here's some code that iterates tries different smoothing parameters until the scipy gives up and raises a RuntimeWarning:
import math
import itertools
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.interpolate import splrep, splev
def is_sorted(a):
return np.all(a[:-1] <= a[1:])
def fit_spline(x, y):
assert(is_sorted(x))
best_n_nans = 2*len(x)
for log_s in itertools.count(0.01, 0.0001):
try:
spl = splrep(x, y, s=math.exp(log_s), k=1)
except RuntimeWarning:
break
n_nans = sum(math.isnan(t) for t in spl[1])
if n_nans < best_n_nans:
best_fit = spl
best_n_nans = n_nans
print(f"Found fit with {n_nans} nans")
if len(spl[0]) <= 5:
break
return best_fit
matplotlib.use("Qt5Agg")
plot = plt.figure()
x = np.array([ 0. , 0. , 0. , 0. , 0. , 0.5, 1.2, 2.3, 2.6,
3. , 3.1, 3.1, 3.1, 5. , 5.1, 5.2, 5.6, 6.3,
6.5, 6.9, 9.6, 9.8, 10. , 10.3, 10.7, 11.3, 11.3,
12.4, 12.4, 12.4, 12.5, 12.5, 13.1, 14.2, 15.4, 15.7,
16.1, 17.7, 18.5, 18.6, 20.5, 21.3, 22.7, 23. , 24.3,
26.3, 26.3, 26.6, 27. , 27.7, 28.7, 31. , 31. , 31.1,
31.1, 31.6, 32.7, 34.4, 34.5, 35.7, 35.7, 38.7, 39.8,
43.9, 43.9, 46.9, 51. , 52.6, 52.6, 53.3, 53.7, 54.5,
57.7, 61.2, 62.5, 65. , 66. , 68.1, 70.8, 74.3, 74.6,
76.3, 83.1, 85.7, 87.5, 90.8, 91. , 92.2, 94.3, 94.6,
101.1, 104.8, 109.6, 113.3, 114.8, 119. , 120.5, 121.1, 126.1,
133.2, 136.2, 137.1, 138.5, 141.4, 144.6, 145.1, 145.3, 150.5,
160.7, 167.4, 168.8, 170.5, 170.5, 173.3, 177. , 181.3, 185.1,
193. , 200.3, 201.8, 204.3, 208.4, 209.9, 223.7, 230.9, 240.3,
243.9, 250.6, 263.1])
y = np.array([ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 6,
0, 0, 0, 9, 0, 0, 0, 0, 0, 0, 9, 1, 0, 3, 0, 9, 2,
10, 0, 14, 19, 7, 10, 18, 22, 31, 12, 29, 16, 20, 2, 38, 26, 28,
26, 25, 17, 35, 20, 29, 39, 51, 29, 27, 43, 38, 40, 47, 41, 39, 54,
44, 53, 65, 50, 59, 51, 62, 69, 74, 66])
spl = fit_spline(x, y)
y2 = splev(x, spl)
plt.plot(x, y, 'o', x, y2)
plot.show()
The result has lots of NaN coefficients:
array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, 35., 38.,
40., 47., 41., 39., 54., 44., 53., 65., 50., 59., 51., 62., 69.,
74., 66., 0., 0.])
Here's a plot of the input data for this example with the partial fit I get:
What can I do to get a reasonable fit? I have tried supplying a weight vector with high weights for the first points but it doesn't do anything. I am not married to scipy if there's another library that does this better.
Ideally I would raise the smoothing parameter until I get a fit with only a few knots, but first I need to get any working fit at all.
BTW, if I supply the knot point the fit works without a problem:
spl = splrep(x, y, w=w, s=math.exp(log_s), k=1, t=[75.0])
gives
However, I'd like to turn this into a repeatable process that doesn't require a human in the loop.