2013-05-01 20 views
5

我有一個一維數據,並與樣條擬合。然後我想找到拐點(忽略鞍點)。現在我正在通過scipy.signal.argrelmin(和argrelmax)搜索splev生成的很多值來搜索它的第一個派生的極值。發現樣條拐點符合1d數據

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() 

但這感覺非常錯誤。我想有可能找到我所需要的而不會產生這麼多的價值。也許像sproot,但適用於第二個派生也許?

+0

肯定的。 ;-) – 2013-06-13 07:25:45

回答

4

The derivative of a B-spline is also a B-spline。因此,您可以首先將樣條擬合到數據中,然後使用導數公式構造導數樣條的係數,最後使用樣條根找出來得到導數樣條的根。這些就是原始曲線的最大值/最小值。

這裏是代碼做到這一點:https://gist.github.com/pv/5504366

係數的相關計算是:

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 

Finding inflection points of a curve via derivative splines

+1

哇,太好了,謝謝。我不知道,區分樣條非常「容易」。由於roots()僅適用於3階,所以如果我想找到極值和拐點,我想我必須使用不同的樣條線,但這對我來說目前不是問題。如果有人對我的問題如何解決您的代碼感興趣,請訪問:https://ideone.com/qKja7X和http://s21.postimg.org/mbc96qcxj/out.png和https://ideone.com/ m757q9 – 2013-05-03 06:53:43