给定一个近似起点,您可以使用一个简单的算法来找到最接近该点的局部最大值。您的拟合代码可能已经这样做了(我不确定您是否成功使用它)。
下面是一些代码,演示了从用户给定的起点进行简单的峰值查找:
#!/usr/bin/env python
from __future__ import division
import numpy as np
from matplotlib import pyplot as plt
# Sample data with two peaks: small one at t=0.4, large one at t=0.8
ts = np.arange(0, 1, 0.01)
xs = np.exp(-((ts-0.4)/0.1)**2) + 2*np.exp(-((ts-0.8)/0.1)**2)
# Say we have an approximate starting point of 0.35
start_point = 0.35
# Nearest index in "ts" to this starting point is...
start_index = np.argmin(np.abs(ts - start_point))
# Find the local maxima in our data by looking for a sign change in
# the first difference
# From http://stackoverflow.com/a/9667121/188535
maxes = (np.diff(np.sign(np.diff(xs))) < 0).nonzero()[0] + 1
# Find which of these peaks is closest to our starting point
index_of_peak = maxes[np.argmin(np.abs(maxes - start_index))]
print "Peak centre at: %.3f" % ts[index_of_peak]
# Quick plot showing the results: blue line is data, green dot is
# starting point, red dot is peak location
plt.plot(ts, xs, '-b')
plt.plot(ts[start_index], xs[start_index], 'og')
plt.plot(ts[index_of_peak], xs[index_of_peak], 'or')
plt.show()
这种方法只有在从您的起点开始完全顺利地上升到顶峰时才有效。如果这需要对噪音更有弹性,我没有使用它,但PyDSTool似乎它可能会有所帮助。这篇SciPy 文章详细介绍了如何使用它来检测嘈杂数据集中的一维峰值。
所以假设此时你已经找到了峰值的中心。现在对于宽度:您可以使用多种方法,但最简单的可能是“半高全宽”(FWHM)。同样,这很简单,因此很脆弱。它会因接近双峰或噪声数据而中断。
FWHM 正是它的名字所暗示的:你会发现峰的宽度是最大值的一半。这是一些执行此操作的代码(它只是从上面继续):
# FWHM...
half_max = xs[index_of_peak]/2
# This finds where in the data we cross over the halfway point to our peak. Note
# that this is global, so we need an extra step to refine these results to find
# the closest crossovers to our peak.
# Same sign-change-in-first-diff technique as above
hm_left_indices = (np.diff(np.sign(np.diff(np.abs(xs[:index_of_peak] - half_max)))) > 0).nonzero()[0] + 1
# Add "index_of_peak" to result because we cut off the left side of the data!
hm_right_indices = (np.diff(np.sign(np.diff(np.abs(xs[index_of_peak:] - half_max)))) > 0).nonzero()[0] + 1 + index_of_peak
# Find closest half-max index to peak
hm_left_index = hm_left_indices[np.argmin(np.abs(hm_left_indices - index_of_peak))]
hm_right_index = hm_right_indices[np.argmin(np.abs(hm_right_indices - index_of_peak))]
# And the width is...
fwhm = ts[hm_right_index] - ts[hm_left_index]
print "Width: %.3f" % fwhm
# Plot to illustrate FWHM: blue line is data, red circle is peak, red line
# shows FWHM
plt.plot(ts, xs, '-b')
plt.plot(ts[index_of_peak], xs[index_of_peak], 'or')
plt.plot(
[ts[hm_left_index], ts[hm_right_index]],
[xs[hm_left_index], xs[hm_right_index]], '-r')
plt.show()
它不一定是半峰全宽——正如一位评论者指出的那样,您可以尝试找出操作员用于峰值检测的正常阈值在哪里,并将其转化为该过程的这一步的算法。
一种更稳健的方法可能是将高斯曲线(或您自己的模型)拟合到以峰值为中心的数据子集 - 例如,从一侧的局部最小值到另一侧的局部最小值 - 并使用其中一个该曲线的参数(例如西格玛)来计算宽度。
我意识到这是很多代码,但我故意避免将索引查找函数分解为“显示我的工作”多一点,当然绘图函数只是为了演示。
希望这至少能给你一个好的起点,让你想出更适合你的特定系列的东西。