在样条拟合的1d数据中找到拐点 [英] finding inflection points in spline fitted 1d data

问题描述

我有一些一维数据，并用样条曲线拟合.然后，我想在其中找到拐点(忽略鞍点).现在，我使用scipy.signal.argrelmin(和argrelmax)搜索由splev生成的许多值，以求其一阶导数的极值.

I have some one dimensional data and fit it with a spline. Then I want to find the inflection points (ignoring saddle points) in it. Now I am searching the extrema of its first derivation by using scipy.signal.argrelmin (and argrelmax) on a lot of values generated by 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()

但这感觉很不对劲.我猜有可能找到我要寻找的东西而不会产生那么多值.像sproot这样的东西，但也许适用于第二推导吗?

But this feels terribly wrong. I guess there is a possibility to find what I am looking for without generating so many values. Something like sproot but applicable to the second derivation perhaps?

推荐答案

The derivative of a B-spline is also a B-spline. You can therefore first fit a spline to your data, then use the derivative formula to construct the coefficients of the derivative spline, and finally use the spline root finding to get the roots of the derivative spline. These are then the maxima/minima of the original curve.

这是执行此操作的代码: https://gist.github.com/pv/5504366

Here is code to do it: 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

这篇关于在样条拟合的1d数据中找到拐点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！