如何使用scipy.curve_fit拟合log(a-x)类型的函数？

Question

我正在尝试拟合一个类似于log(y)=a*log(b-x)+c的函数，其中a、b和c是需要拟合的参数。相关的代码如下：

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

def logfunc(T, a, b, c):
    v=(a*np.log(b-T))+c
    return v

popt, pcov=curve_fit(logfunc, T, np.log(Energy), check_finite=False, bounds=([0.1, 1.8, 0.1], [1.0, 2.6, 1.0]))

plt.plot(T, logfunc(T, *popt))
plt.show

这里的T和Energy是一些生成的数据（我用它来绘制其他东西，所以数据应该没问题）。T在0.3到3.2之间。我很确定问题出在b=T这个点上，因为我一直收到错误信息ValueError: Residuals are not finite in the initial point，但我不知道如何解决。

- Ilin Karagjozov

---

永远不要假设数据是正确的。即使在其他地方它是OK的，也可能在这里导致问题。最好直接在示例中创建人工数据。这使得问题对他人来说可以验证（而且作为额外的好处，找出哪个特定的数据组合导致问题可能会带你找到实际的解决方案）。 - MB-F

需要注意的是：如果 T 在 0.3 和 3.2 之间，当 b 被限制在 1.8 和 2.8 之间时，你认为 log(b-T) 的结果会是什么？（同时我不确定初始值是否都为 1，无论边界如何。） - MB-F

@kazemakase 这是一个重要的观点。无论如何，设置 p0 可能是一个好主意。此外，如果 min(b) < max(T)，就会遇到问题。因此，应该检查这一点。最后的问题是：为什么不适合指数版本呢？这样就可以避免 log( negativNumber ) 问题。 - mikuszefski

@kazemakase 謝謝您的回應。我比較新於在這裡發問，因此對禮儀不太熟悉，所以我會記住您關於數據的建議。 - Ilin Karagjozov

@IlinKaragjozov 这部分是礼节，但更多的是为了你自己的利益。你越能让别人轻松地复制和修复你的问题，就越有可能得到答案。如果别人可以简单地复制、粘贴和运行代码，那么人们就更有可能深入研究问题。 - MB-F

Answer 1

您可能会发现 lmfit 包 (http://lmfit.github.io/lmfit-py/) 对于这类问题非常有用。它提供了一种更高级的曲线拟合方法，以及比 scipy.optimize 包或 curve_fit() 函数更好的参数和模型抽象化。

对于这个问题，lmfit 的两个重要特点是：

1. 能够设置变量的边界。 curve_fit() 也可以做到这一点，但只能通过使用有序列表来设置最小/最大边界。而在 lmfit 中，边界属于 Parameter 对象。
2. 有一种明确设置处理 NaN 值策略的方法，这可能会对您的拟合造成问题。

使用 lmfit，您的脚本将编写为：

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model

def logfunc(T, a, b, c):
    return (a*np.log(b-T))+c

log_model = Model(logfunc, nan_policy='raise')  # raise error on NaNs
params = log_model.make_params(a=0.5, b=2.0, c=0.5) # initial values
params['b'].min = 1.8  # set min/max values
params['b'].max = 2.6 
params['c'].min = 0.1  # and so forth 

result = log_model.fit(np.log(Energy), params, T=T)

print(result.fit_report())

plt.plot(T, Energy, 'bo', label='data')
plt.plot(T, np.exp(result.best_fit), 'r--', label='fit')
plt.legend()
plt.xlabel('T')
plt.ylabel('Energy')
plt.gca().set_yscale('log', basey=10)
plt.show()

这个脚本比你的起始脚本稍微冗长一些，因为它提供了一个带标签的图表，而且使用参数对象而不是标量可以提供更多的灵活性和清晰度。

对于您的拟合，您可能考虑将nan_policy设置为'omit'，这样会在出现NaN时省略它们——虽然这从来不是一个好主意，但有时有助于帮助您找到log(b-T)有效的位置。您还可以修改模型函数以执行类似的操作

def logfunc(T, a, b, c):
    arg = b - T
    arg[np.where(arg < 1.e-16)] = 1.e-16
    return a*np.log(arg) + c

为明确防止NaN的一个明显原因。

- M Newville

Answer 2

残差在初始点不是有限的

这意味着初始点不好，其中一些对数是无限或未定义的。您需要一个更好的初始点。

由于模型的特性，b必须大于T中的任何点。您目前设置的b边界并不能保证这一点。请加强它们。

当您不提供p0参数时，SciPy会在提供的边界内猜测。因此，如果边界保证有限性，则不会发生错误。尽管如此，通常最好自己指定p0，因为您对问题具有更好的先验了解，而SciPy则没有。

带有调整边界的工作示例：

popt, pcov=curve_fit(logfunc, np.linspace(0.3, 3.2, 6), [8, 7, 6, 5, 4, 3], bounds=([0.1, 3.2, 0.1], [1.0, 3.6, 1.0]))

- user6655984