1

我有以下 Python (v2.7.14) 代码，它使用来自 SciPy (v1.0.1) 的 curve_fit 来查找指数衰减函数的参数。大多数时候，我得到合理的结果。有时，我会得到一些完全超出我预期范围的结果，即使找到的参数在与原始图表进行绘制时看起来很好。

首先，我对指数衰减公式的理解来自https://en.wikipedia.org/wiki/Exponential_decay，我将其翻译成 Python 如下：

y = a * numpy.exp(-b * x) + c

在哪里：

  • a是数据的初始值
  • b是衰减率，它是信号从初始值到 1/e 时的倒数
  • c是一个偏移量，因为我正在处理我的数据中永远不会达到零的非负值
  • x是当前时间

该脚本考虑到正在拟合非负数据并适当地抵消初始猜测。但即使没有猜测，没有偏移，使用最大/最小（而不是第一个/最后一个值）和我尝试过的其他随机事物，我似乎也无法让 curve_fit 在麻烦的数据集上产生合理的值。

我的假设是，麻烦的数据集没有足够的曲线可以在不超出数据范围的情况下进行拟合。我查看了 curve_fit 的 bounds 参数，并认为这可能是一个合理的选择。我不确定什么是计算的好下限和上限，或者它是否真的是我正在寻找的选项。

这是代码。注释掉的代码是我尝试过的东西。

#!/usr/local/bin/python

import numpy as numpy
from scipy.optimize import curve_fit
import matplotlib.pyplot as pyplot

def exponential_decay(x, a, b, c):
    return a * numpy.exp(-b * x) + c

def fit_exponential(decay_data, time_data, decay_time):
    # The start of the curve is offset by the last point, so subtract
    guess_a = decay_data[0] - decay_data[-1]
    #guess_a = max(decay_data) - min(decay_data)

    # The time that it takes for the signal to reach 1/e becomes guess_b
    guess_b = 1/decay_time

    # Since this is non-negative data, above 0, we use the last data point as the baseline (c)
    guess_c = decay_data[-1]
    #guess_c = min(decay_data)

    guess=[guess_a, guess_b, guess_c]
    print "guess: {0}".format(guess)

    #popt, pcov = curve_fit(exponential_decay, time_data, decay_data, maxfev=20000)
    popt, pcov = curve_fit(exponential_decay, time_data, decay_data, p0=guess, maxfev=20000)

    #bound_lower = [0.05, 0.05, 0.05]
    #bound_upper = [decay_data[0]*2, guess_b * 10, decay_data[-1]]
    #print "bound_lower: {0}".format(bound_lower)
    #print "bound_upper: {0}".format(bound_upper)
    #popt, pcov = curve_fit(exponential_decay, time_data, decay_data, p0=guess, bounds=[bound_lower, bound_upper], maxfev=20000)

    a, b, c = popt

    print "a: {0}".format(a)
    print "b: {0}".format(b)
    print "c: {0}".format(c)

    plot_fit = exponential_decay(time_data, a, b, c)

    pyplot.plot(time_data, decay_data, 'g', label='Data')
    pyplot.plot(time_data, plot_fit, 'r', label='Fit')
    pyplot.legend()
    pyplot.show()

print "Gives reasonable results"
time_data = numpy.array([0.0,0.040000000000000036,0.08100000000000018,0.12200000000000011,0.16200000000000014,0.20300000000000007,0.2430000000000001,0.28400000000000003,0.32400000000000007,0.365,0.405,0.44599999999999995,0.486,0.5269999999999999,0.567,0.6079999999999999,0.6490000000000002,0.6889999999999998,0.7300000000000002,0.7700000000000002,0.8110000000000002,0.8510000000000002,0.8920000000000001,0.9320000000000002,0.9730000000000001])
decay_data = numpy.array([1.342146870531986,1.405586070225509,1.3439802492549762,1.3567811728250267,1.2666276377825874,1.1686375326985337,1.216119360088685,1.2022841507836042,1.1926979408026064,1.1544395213303447,1.1904416926531907,1.1054720201415882,1.112100683833435,1.0811434035632939,1.1221671794680403,1.0673295063196415,1.0036146509494743,0.9984005680821595,1.0134498134883763,0.9996920772051201,0.929782730581616,0.9646581154122312,0.9290690593684447,0.8907360533169936,0.9121560047238627])
fit_exponential(decay_data, time_data, 0.567)

print

print "Gives results that are way outside my expectations"
time_data = numpy.array([0.0,0.040000000000000036,0.08099999999999996,0.121,0.16199999999999992,0.20199999999999996,0.24300000000000033,0.28300000000000036,0.32399999999999984,0.3650000000000002,0.40500000000000025,0.44599999999999973,0.48599999999999977,0.5270000000000001,0.5670000000000002,0.6079999999999997,0.6479999999999997,0.6890000000000001,0.7290000000000001,0.7700000000000005,0.8100000000000005,0.851,0.8920000000000003,0.9320000000000004,0.9729999999999999,1.013,1.0540000000000003])
decay_data = numpy.array([1.4401611921948776,1.3720688158534153,1.3793465463227048,1.2939909686762128,1.3376345321949346,1.3352710161631154,1.3413634841956348,1.248705138603995,1.2914294791901497,1.2581763134585313,1.246975264018646,1.2006447776495062,1.188232179689515,1.1032789127515186,1.163294324147017,1.1686263160765304,1.1434009568472243,1.0511578409946472,1.0814520440570896,1.1035953824496334,1.0626893599266163,1.0645580326776076,0.994855722989818,0.9959891485338087,0.9394584009825916,0.949504060086646,0.9278639431146273])
fit_exponential(decay_data, time_data, 0.6890000000000001)

这是文本输出：

Gives reasonable results
guess: [0.4299908658081232, 1.7636684303350971, 0.9121560047238627]
a: 1.10498934435
b: 0.583046565885
c: 0.274503681044

Gives results that are way outside my expectations
guess: [0.5122972490802503, 1.4513788098693758, 0.9278639431146273]
a: 742.824622191
b: 0.000606308344957
c: -741.41398516

最值得注意的是，对于第二组结果，a的值非常高，c的值在负值范围内同样低，b是一个非常小的十进制数。

这是第一个数据集的图表，它给出了合理的结果。 这是第一个数据集的图表，它给出了合理的结果。

这是第二个数据集的图表，它没有给出好的结果。 这是第二个数据集的图表，它没有给出好的结果。

请注意，图形本身绘制正确，尽管该线实际上并没有很好的曲线。

我的问题：

  • 我用curve_fit实现的指数衰减算法是否正确？
  • 我的初始猜测参数是否足够好？
  • bounds 参数是这个问题的正确解决方案吗？如果是这样，确定下限和上限的好方法是什么？
  • 我在这里错过了什么吗？

再次谢谢你！

4

1 回答 1

0

当您说第二次拟合给出的结果“超出”您的预期，并且尽管第二张图“绘制正确”时，这条线并没有真正“具有良好的曲线拟合”，您在理解什么是正确的轨道上继续。我认为你只是错过了一块拼图。

第二张图与看起来是线性的曲线非常吻合。这可能意味着您的数据实际上没有足够的变化（嗯，可能低于噪声水平）来检测它是指数衰减。

我敢打赌，如果您不仅打印出最佳拟合值，还打印出变量的不确定性和相关性，您会发现不确定性很大，并且一些相关性非常接近 1。这可能意味着考虑考虑不确定性（并且测量总是有不确定性），结果可能实际上符合您的期望。这也可能告诉您，您拥有的数据不能很好地支持指数衰减。

您也可以为这些数据尝试其他模型（想到“线性”；））并比较拟合优度统计数据，例如卡方和 Akaike 信息标准。

scipy.curve_fit确实返回协方差矩阵 -pcov您在示例中未使用的。不幸的是，scipy.curve_fit不会将这些值转换为不确定性和相关值，并且根本不会尝试返回任何拟合优度统计信息。

要完全解释对数据的任何拟合，您不仅需要最佳拟合值，还需要对可变参数的不确定性进行估计。您需要拟合优度统计数据来确定拟合是否良好，或者至少一种拟合是否优于另一种。

于 2018-07-17T01:33:35.463 回答