znajdowanie punktów przegięcia w dopasowanych splajtach danych 1d

Question

Mam dane jednowymiarowe i dopasowuję je za pomocą splajnu. Następnie chcę znaleźć w nim punkty przegięcia (ignorując punkty siodełka). Teraz przeszukuję extrema jego pierwszej derywacji za pomocą scipy.signal.argrelmin (i argrelmax) na wielu wartościach generowanych przez 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() 

Ale to strasznie źle. Sądzę, że istnieje możliwość znalezienia tego, czego szukam, bez generowania tak wielu wartości. Coś jak sproot, ale może być zastosowane do drugiej derywacji?

Answer 1

The derivative of a B-spline is also a B-spline. Możesz zatem najpierw dopasować splajn do swoich danych, a następnie użyć formuły pochodnej do skonstruowania współczynników pochodnej splajnu, a na końcu użyć znalezienia splajnu, aby uzyskać korzenie pochodnej splajnu. Są to wówczas maksima/minima pierwotnej krzywej.

Oto kod, aby to zrobić: https://gist.github.com/pv/5504366

Odpowiednie obliczenia współczynników jest:

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