python - 使用边界时 scipy.optimize.curvefit 失败

Question

我正在尝试使用scipy.optimize.curvefit将一组数据与函数相匹配(请参见下面的示例)， 但是当我使用边界( documentation )时，拟合失败，我只是得到 初始猜测参数作为输出。 一旦我替换 -np.inf ad np.inf 作为第二个参数的边界 (函数中的dt)，拟合有效。 我做错了什么？

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt

#Generate data
crc=np.array([-1.4e-14, 7.3e-14, 1.9e-13, 3.9e-13, 6.e-13, 8.0e-13, 9.2e-13, 9.9e-13, 
              1.e-12, 1.e-12, 1.e-12, 1.0e-12, 1.1e-12, 1.1e-12, 1.1e-12, 1.0e-12, 1.1e-12])
time=np.array([0., 368., 648., 960., 1520.,1864., 2248., 2655., 3031., 
                3384., 3688., 4048., 4680., 5343., 6055.,  6928.,  8120.])

#Define the function for the fit
def testcurve(x, Dp, dt):
    k = -Dp*(x+dt)*2e11
    curve = 1e-12 * (1+2*(-np.exp(k) + np.exp(4*k) - np.exp(9*k) + np.exp(16*k)))
    curve[0]= 0
    return curve

#Set fit bounds 
dtmax=time[2]
param_bounds = ((-np.inf, -dtmax),(np.inf, dtmax))

#Perform fit
(par, par_cov) = opt.curve_fit(testcurve, time, crc, p0 = (5e-15, 0), bounds = param_bounds)

#Print and plot output
print(par)       
plt.plot(time, crc, 'o')
plt.plot(time, testcurve(time, par[0], par[1]), 'r-')
plt.show()

Answer 1

我今天在不同的拟合问题中遇到了相同的行为。经过一番网上搜索，我发现这个链接很有帮助:Why does scipy.optimize.curve_fit not fit to the data?

简短的回答是:在数值拟合中使用极小(或大)的数字并不稳健，缩放它们会带来更好的拟合。

<小时/>

就您而言，crc 和 Dp 都是非常小的数字，可以按比例放大。您可以使用比例因子，在一定范围内，拟合看起来相当稳健。完整示例:

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt

#Generate data
crc=np.array([-1.4e-14, 7.3e-14, 1.9e-13, 3.9e-13, 6.e-13, 8.0e-13, 9.2e-13, 9.9e-13, 
              1.e-12, 1.e-12, 1.e-12, 1.0e-12, 1.1e-12, 1.1e-12, 1.1e-12, 1.0e-12, 1.1e-12])
time=np.array([0., 368., 648., 960., 1520.,1864., 2248., 2655., 3031., 
                3384., 3688., 4048., 4680., 5343., 6055.,  6928.,  8120.])

# add scale factors to the data as well as the fitting parameter
scale_factor_1 = 1e12 # 1./np.mean(crc) also works if you don't want to set the scale factor manually
scale_factor_2 = 1./2e11

#Define the function for the fit
def testcurve(x, Dp, dt):
    k = -Dp*(x+dt)*2e11 * scale_factor_2
    curve = 1e-12 * (1+2*(-np.exp(k) + np.exp(4*k) - np.exp(9*k) + np.exp(16*k))) * scale_factor_1
    curve[0]= 0
    return curve

#Set fit bounds 
dtmax=time[2]
param_bounds = ((-np.inf, -dtmax),(np.inf, dtmax))

#Perform fit
(par, par_cov) = opt.curve_fit(testcurve, time, crc*scale_factor_1, p0 = (5e-15/scale_factor_2, 0), bounds = param_bounds)

#Print and plot output
print(par[0]*scale_factor_2, par[1])

plt.plot(time, crc*scale_factor_1, 'o')
plt.plot(time, testcurve(time, par[0], par[1]), 'r-')
plt.show()

拟合结果:[6.273102923176595e-15, -21.12202697564494]，它提供了合理的拟合，并且非常接近没有任何边界的结果:[6.27312512e-15, - 2.11307470e+01]