Finding the smoothness of a spline using scipy

Question

Consider the following Example:

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy import interpolate
xs = np.linspace(1,10,500)
ys = [0.92 * x ** 2.3 + 0.0132 * x ** 4 + 0.0743 * (x - 9) ** 3 - 4 * (x -3) ** 2 + 80 * math.sin(math.sin(x)) + 10 * math.sin(x*5) + 1.2* np.random.normal(-4,4,1) for x in xs]
ys[200] = ys[200] + 130
ys[201] = ys[201] + 135
ys[202] = ys[202] + 129
ys[203] = ys[203] + 128
ys[204] = ys[204] + 131
ys[205] = ys[205] + 130
ys[206] = ys[206] + 129
ys[207] = ys[207] + 129
ys[208] = ys[208] + 128
ys[209] = ys[209] + 130

If I plot xs and ys at this point, it produces a nice graph:

Now I am using scipy.interpolate.splrep to fit a spline curve to this data. I have used two different splines to fit two different segments of the data:

tck = interpolate.splrep(xs[0:199], ys[0:199], s = 1000)
ynew2 = interpolate.splev(xs[0:199], tck, der = 0)

and :

tck = interpolate.splrep(xs[210:500], ys[210:500], s = 9000)
ynew3 = interpolate.splev(xs[210:500], tck, der = 0)

Then we have :

Now I want to programmatically detect the quality of the fit. The fit should neither be too straight - i.e. preserve features, nor should it "overdetect" the noisy variations as features.

I plan to use a peak counter fed to an ANN.

However, at this point, my question is:

Does scipy/numpy have a built in function where i can feed in the output of splrep and it will compute the minima or maxima and the density of maxima/minima at any particular interval?

Note:
I am aware of the R**2 value, I am looking to find another measure to detect preservation of features.

Unknown Unknown · Accepted Answer · 2017-12-18T23:06:10

SciPy does not have a method for finding the critical points of a cubic spline. The closest we have sproot which find the roots of a cubic spline. In order for this to be useful here, we must fit splines of order 4, so that the derivative is a cubic spline. This is what I do below

from scipy.interpolate import splrep, splev, splder, sproot

tck1 = splrep(xs[0:199], ys[0:199], k=4, s=1000)
tck2 = splrep(xs[210:500], ys[210:500], k=4, s=9000)
roots1 = sproot(splder(tck1), 1000)     # 1000 is an upper bound for the number of roots
roots2 = sproot(splder(tck2), 1000)

x1 = np.linspace(xs[0], xs[198], 1000)     # plot both splines
plt.plot(x1, splev(x1, tck1))
x2 = np.linspace(xs[210], xs[499], 1000)
plt.plot(x2, splev(x2, tck2))             

plt.plot(roots1, splev(roots1, tck1), 'ro')        # plot their max/min points
plt.plot(roots2, splev(roots2, tck2), 'ro') 
plt.show()

The difference is obvious.

You can also find the number of roots in any specific interval, such as [3, 4]:

np.where((3 <= roots1) & (roots1 <= 4))[0].size    # 29

or equivalently, np.sum((3 <= roots1) & (roots1 <= 4))

Finding the smoothness of a spline using scipy

1 Answers