我正在尝试计算两种形式的指数对某些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的奇怪结果';s(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]
所以现在问题变为:我如何才能自动为第二个函数提供更合理的参数?
Zunzun.com使用差分进化遗传算法(DE)找到初始参数估计值,然后将其传递给scipy中的Levenberg-Marquardt求解器。DE本身并不是一个全局优化器,而是一个"初始参数猜测器"。
你可以在网站的任何网页底部找到zunzun.com安装程序的BSD许可Python源代码链接——它有很多全面的例子——所以现在不需要自己编写代码。如果你有任何问题,请告诉我,我会尽我所能提供帮助。
詹姆斯·菲利普斯zunzun@zunzun.com
请注意,在zunzun估计的a=0
和您的第一个模型中。所以他们只是在估计一个常数。因此,第一种情况下的b
和第二种情况下b
和c
是不相关的,并且没有被识别。
Zunzun还使用差分进化作为全局求解器,这是我上次看到的。Scipy现在使用basinhoping作为全局优化器,看起来非常好,在可能出现局部极小值的情况下值得一试。
我的"廉价"方式,因为参数在你的例子中没有很大的范围:尝试随机起始值
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为尊尊,cf为曲线拟合)
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
和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
第二个函数与第一个函数没有什么不同。CCD_ 13和CCD_。因此,使用导数信息的优化器也会遇到问题,因为如果我没有错的话,雅可比在某些方向上是奇异的。