Trouver la douceur d'une spline en utilisant scipy

Question

Considérez l'exemple suivant:

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

Si je trace xs et ys à ce stade, cela produit un joli graphique:

Maintenant j'utilise scipy.interpolate.splrep pour ajuster une courbe spline à ces données. J'ai utilisé deux splines différentes pour ajuster deux segments différents des données:

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

et :

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

Ensuite nous avons :

Maintenant je veux détecter de manière programmative la qualité de l'ajustement. L'ajustement ne doit pas être trop droit - c'est-à-dire préserver les caractéristiques, ni "surestimer" les variations bruyantes comme des caractéristiques.

Je prévois d'utiliser un compteur de pics alimenté à un ANN.

Cependant, à ce stade, ma question est:

• Est-ce que scipy/numpy a une fonction intégrée où je peux fournir la sortie de splrep et elle calculera les minima ou maxima et la densité des maxima/minima à un intervalle particulier?

Note:
Je suis conscient de la valeur R**2, je cherche à trouver une autre mesure pour détecter la préservation des caractéristiques.

Answer 1

SciPy n'a pas de méthode pour trouver les points critiques d'une spline cubique. Le plus proche que nous avons est sproot qui trouve les racines d'une spline cubique. Pour que cela soit utile ici, nous devons ajuster des splines d'ordre 4, de sorte que la dérivée soit une spline cubique. C'est ce que je fais ci-dessous

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 est une borne supérieure pour le nombre de racines
roots2 = sproot(splder(tck2), 1000)

x1 = np.linspace(xs[0], xs[198], 1000)     # tracer les deux 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')        # tracer leurs points max/min
plt.plot(roots2, splev(roots2, tck2), 'ro') 
plt.show()

La différence est évidente.

Vous pouvez également trouver le nombre de racines dans n'importe quel intervalle spécifique, comme [3, 4]:

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

ou de manière équivalente, np.sum((3 <= roots1) & (roots1 <= 4))