2

我正在尝试计算两种形式的指数对某些x, y数据的最佳拟合（数据文件可以从这里下载）

这是代码：

from scipy.optimize import curve_fit
import numpy as np

# Get x,y data
data = np.loadtxt('data.txt', unpack=True)
xdata, ydata = data[0], data[1]

# Define first exponential function
def func(x, a, b, c):
    return a * np.exp(b * x) + c

# Get parameters estimate
popt, pcov = curve_fit(func, xdata, ydata)

print popt

# Define second exponential function (one more parameter)
def func2(x, a, b, c, d):
    return a * np.exp(b * x + c) + d

# Get parameters estimate
popt2, pcov2 = curve_fit(func2, xdata, ydata)

print popt2

第一个指数给出的值与 zunzun.com（此处为 PDF）完全相同popt：

[  7.67760545e-15   1.52175476e+00   2.15705939e-02]

但第二个给出的值显然是错误的popt2：

[ -1.26136676e+02  -8.13233297e-01  -6.66772692e+01   3.63133641e-02]

这是同一第二个功能的 zunzun.com 值（此处为 PDF）：

a = 6.2426224704624871E-15
b = 1.5217697532005228E+00
c = 2.0660424037614489E-01
d = 2.1570805929514186E-02

我尝试按照此处推荐的方式制作列表数组使用 python (scipy) 曲线拟合的奇怪结果，但这并没有帮助。我在这里做错了什么？

添加 1

我猜这个问题与我提供函数时缺少初始值有关（如此处所述：高斯拟合与 python 中的 scipy.optimize.curve_fit 结果错误）

如果我像这样将第一个指数的估计值提供给第二个指数（使新参数d最初为零）：

popt2, pcov2 = curve_fit(func2, xdata, ydata, p0 = [popt[0], popt[1], popt[2], 0])

与 zunzun.com 相比，我得到的结果非常合理但仍然错误：

[  1.22560853e-14   1.52176160e+00  -4.67859961e-01   2.15706930e-02]

所以现在问题变成了：如何自动为我的第二个函数提供更合理的参数？

4

2 回答 2

3

Zunzun.com 使用差分进化遗传算法 (DE) 来查找初始参数估计，然后将其传递给 scipy 中的 Levenberg-Marquardt 求解器。DE 本身实际上并没有用作全局优化器，而是用作“初始参数猜测器”。

您可以在该站点的任何网页底部找到 zunzun.com fitter 的 BSD 许可 Python 源代码的链接 - 它有许多全面的示例 - 因此无需立即自己编写代码。如果您有任何问题，请告诉我，我会尽力提供帮助。

詹姆斯·菲利普斯 zunzun@zunzun.com

于 2013-07-09T14:12:07.530 回答
1

请注意，a=0在 zunzun 的估计和您的第一个模型中。所以他们只是在估计一个常数。因此，b在第一种情况下和b在c第二种情况下是不相关的并且没有被识别。

Zunzun 还使用微分进化作为全局求解器，这是我上次看的时候。Scipy 现在将盆地跳跃作为全局优化器，看起来相当不错，在可能存在局部最小值的情况下值得一试。

我的“便宜”方式，因为参数在您的示例中没有很大的范围：尝试随机起始值

np.random.seed(1)
err_last = 20
best = None

for i in range(10):
    start = np.random.uniform(-10, 10, size=4)
    # Get parameters estimate
    try:
        popt2, pcov2 = curve_fit(func2, xdata, ydata, p0=start)
    except RuntimeError:
        continue
    err = ((ydata - func2(xdata, *popt2))**2).sum()
    if err < err_last:
        err_last = err
        print err
        best = popt2


za = 6.2426224704624871E-15
zb = 1.5217697532005228E+00
zc = 2.0660424037614489E-01
zd = 2.1570805929514186E-02

zz = np.array([za,zb,zc,zd])
print 'zz', zz
print 'cf', best

print 'zz', ((ydata - func2(xdata, *zz))**2).sum()
print 'cf', err_last

最后一部分打印（zz 是zunzun，cf 是curve_fit）

zz [  6.24262247e-15   1.52176975e+00   2.06604240e-01   2.15708059e-02]
cf [  1.24791299e-16   1.52176944e+00   4.11911831e+00   2.15708019e-02]
zz 9.52135153898
cf 9.52135153904

b和的参数与 Zunzun 不同c，但残差平方和相同。

添加

a * np.exp(b * x + c) + d = np.exp(b * x + (c + np.log(a))) + d

或者

a * np.exp(b * x + c) + d = (a * np.exp(c)) * np.exp(b * x) + d

第二个函数与第一个函数并没有什么不同。a并且c没有单独标识。所以使用导数信息的优化器也会有问题，因为雅可比在某些方向上是奇异的，如果我没看错的话。

于 2013-07-09T02:12:42.637 回答