trovare i punti di flesso nella spline montata 1d dati

Question

Ho alcuni dati unidimensionali e li inserisco con una spline. Quindi voglio trovare i punti di flesso (ignorando i punti della sella) in esso. Ora sto cercando gli estremi della sua prima derivazione usando scipy.signal.argrelmin (e argrelmax) su molti valori generati da splev.

import scipy.interpolate 
import scipy.optimize 
import scipy.signal 
import numpy as np 
import matplotlib.pyplot as plt 
import operator 

y = [-1, 5, 6, 4, 2, 5, 8, 5, 1] 
x = np.arange(0, len(y)) 
tck = scipy.interpolate.splrep(x, y, s=0) 

print 'roots', scipy.interpolate.sproot(tck) 
# output: 
# [0.11381478] 

xnew = np.arange(0, len(y), 0.01) 
ynew = scipy.interpolate.splev(xnew, tck, der=0) 

ynew_deriv = scipy.interpolate.splev(xnew, tck, der=1) 

min_idxs = scipy.signal.argrelmin(ynew_deriv) 
max_idxs = scipy.signal.argrelmax(ynew_deriv) 
mins = zip(xnew[min_idxs].tolist(), ynew_deriv[min_idxs].tolist()) 
maxs = zip(xnew[max_idxs].tolist(), ynew_deriv[max_idxs].tolist()) 
inflection_points = sorted(mins + maxs, key=operator.itemgetter(0)) 

print 'inflection_points', inflection_points 
# output: 
# [(3.13, -2.9822449358974357), 
# (5.03, 4.3817785256410255) 
# (7.13, -4.867132628205128)] 

plt.legend(['data','Cubic Spline', '1st deriv']) 
plt.plot(x, y, 'o', 
     xnew, ynew, '-', 
     xnew, ynew_deriv, '-') 
plt.show() 

Ma questo si sente terribilmente sbagliato. Immagino che ci sia la possibilità di trovare quello che sto cercando senza generare così tanti valori. Qualcosa come il germoglio ma applicabile alla seconda derivazione forse?

Answer 1

derivative of a B-spline is also a B-spline. È quindi possibile adattare prima una spline ai dati, quindi utilizzare la formula derivativa per costruire i coefficienti della spline derivata e infine utilizzare la ricerca della radice spline per ottenere le radici della spline derivata. Questi sono quindi i massimi/minimi della curva originale.

Ecco il codice per farlo: https://gist.github.com/pv/5504366

Il calcolo dei coefficienti pertinente è:

t, c, k = scipys_spline_representation 
# Compute the denominator in the differentiation formula. 
dt = t[k+1:-1] - t[1:-k-1] 
# Compute the new coefficients 
d = (c[1:-1-k] - c[:-2-k]) * k/dt 
# Adjust knots 
t2 = t[1:-1] 
# Pad coefficient array to same size as knots (FITPACK convention) 
d = np.r_[d, [0]*k] 
# Done, a new spline 
new_spline_repr = t2, d, k-1