我试图计算两种指数形式对某些x, y数据的最佳拟合(数据文件可以从here下载)
代码如下:

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

对于popt,第一个指数给出的值与zunzun.com(PDF here)完全相同:
[  7.67760545e-15   1.52175476e+00   2.15705939e-02]

但是第二个给出的值显然是错误的:
[ -1.26136676e+02  -8.13233297e-01  -6.66772692e+01   3.63133641e-02]

这是同一个第二个函数的zunzun.com值(PDF here):
a = 6.2426224704624871E-15
b = 1.5217697532005228E+00
c = 2.0660424037614489E-01
d = 2.1570805929514186E-02

我试着在这里创建列表数组Strange result with python's (scipy) curve fitting,但这没有帮助。我在这里做错什么了?
加1
我猜这个问题与我输入函数的初始值不足有关(如这里所述:gaussian fit with scipy.optimize.curve_fit in python with wrong results
如果我将第一个指数的估计值输入到第二个指数中(使新参数popt2初始值为零):
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]

所以现在的问题变成了:我怎样才能自动给我的第二个函数提供更合理的参数?

最佳答案

请注意,在zunzun的估算和您的第一个模型中a=0。所以他们只是估计一个常数。因此,第一种情况下的b和第二种情况下的bc是不相关的,没有被识别。
在我最后一次看的时候,Zunzun也将差分进化作为一个全局解算器。Scipy现在有了basinhopping作为全局优化器,看起来非常不错,在可能出现局部极小值的情况下值得一试。
我的“廉价”方法,因为在您的示例中,参数没有很大的范围:尝试随机起始值

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

与尊尊尊不同的参数有bc,但残差平方和相同。
添加
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
第二个函数和第一个函数没有什么区别。ac不单独标识。所以使用导数信息的优化器也会有问题,因为雅可比在某些方向上是奇异的,如果我看得对的话。

10-07 14:54