在对指数关系进行curve_fit拟合时，某个参数的不确定性非常高
用同一个函数拟合的不同图形
我尝试用同一个函数拟合我的图表，本应呈现负指数型，但对L∞ 的预测却非常差。如何改进曲线拟合？当我绘制L∞ 的误差棒时，预测的不确定性竟然非常大，而其他参数的拟合要好得多。
代码我用的：
def twentyconc_model(x, k, Linfty):
return Linfty - Linfty*exp(-k*x)
popt20 = []
pcov20 = []
for i in range(0,5):
popt, pcov = curve_fit(twentyconc_model, twentyconc_data[i][:,0], twentyconc_data[i][:,1], p0=[0, mean(twentyconc_data[i][-10,1])], maxfev=10000)
popt20.append(popt)
pcov20.append(pcov)
print(popt20,"\n\n",pcov20)
for i in range(0,5):
twentyk = popt20[i][0] # estimated k at different temperatures @ 20 vol%
twentyLinfty = popt20[i][1] # estimated Linfty at different temperatures @ 20 vol%
twentyconc_time_model = np.linspace(min(twentyconc_data[i][:,0]), max(twentyconc_data[i][:,0]), 100)
twentyconc_conduct_model = twentyconc_model(twentyconc_time_model, twentyk, twentyLinfty)
plt.plot(twentyconc_time_model, twentyconc_conduct_model, label = label20[i])
#plot(time_shifted20, conductance20, label = label20[i], linestyle=' ', marker='o')
#plt.legend()
#plt.show()
legend()
show()
解决方案
这里的问题并不是curve_fit在数值上失败，而是从当前数据中很难辨识出L∞，这与这个针对两个强耦合参数的模型类似；该模型的k 控制曲线升高的速度，L∞ 控制渐近平台。由于这里提供的数据并未足够接近平台，许多k 和L∞ 的组合会产生几乎相同的曲线。若存在较大的k 搭配较小的L∞，或较小的k 搭配较大的L∞，这两种情况都能同样好地拟合观测区域，因此导致协方差矩阵接近奇异，L∞ 的不确定性会大幅增加。
本质上，一些需要改进的地方是把对L∞ 的初始猜测从 y[-10,1] 改为 np.mean(y[-10:,1])，并改用最近10次测量的平均值来估算平台值，以获得更现实的渐近值，从而使收敛更稳定。以 k = 0 开始相当于一条平直线，与预期的指数行为相距甚远，因此一个合理的正向初始猜测有助于更快找到正确解（例如；0.02）。参数界限（k > 0, L∞ > 0）将防止优化器在物理上毫无意义的参数值上探索，从而降低不稳定性。如果你使用鲁棒拟合（loss='soft_l1'），它会减少离群值的影响，防止少量异常点对拟合参数产生过大影响。
但除了这些，对于提供的代码，如果你按如下方式修改：
def twentyconc_model(t, k, Linfty):
return Linfty * (1 - np.exp(-k * t))
fit_results = []
for i in range(len(twentyconc_data)):
x = twentyconc_data[i][:, 0]
y = twentyconc_data[i][:, 1]
# Better initial guesses
Linfty_guess = np.mean(y[-10:])
k_guess = 0.02
popt, pcov = curve_fit(twentyconc_model,x,y,p0=[k_guess, Linfty_guess],
bounds=(
[0, 0], # k > 0, Linfty > 0
[np.inf, np.inf]
),maxfev=50000)
k_fit, Linfty_fit = popt
# Standard errors
perr = np.sqrt(np.diag(pcov))
k_err = perr[0]
Linfty_err = perr[1]
# Correlation matrix
corr = pcov / np.outer(perr, perr)
fit_results.append({"k": k_fit,"k_err": k_err,"Linfty": Linfty_fit,"Linfty_err": Linfty_err,"corr": corr})
# Plot fit
x_fit = np.linspace(x.min(), x.max(), 500)
y_fit = twentyconc_model(x_fit, *popt)
plt.scatter(x, y, s=15)
plt.plot(x_fit, y_fit, label=label20[i])
并在 fit_results 处查看数据：
for i, res in enumerate(fit_results):
print(f"\ndataset {i}")
print(f"k= {res['k']:.6f} +- {res['k_err']:.6f}")
print(f"Linfty= {res['Linfty']:.3f} +- {res['Linfty_err']:.3f}")
print("Correlation matrix")
print(res['corr'])
如果k 与L∞ 的相关性仍然很高（即类似于：[[ 1.00 -0.998] [-0.998 1.00]]），根本问题在于数据不足以约束L∞，在这种情况下获取更多接近平台的测量值或固定L∞ 将比再调试优化器更有效。
备选方案1
参数不稳定的主要原因是在三张图中存在那些极端离群值。如果在原始数据集中排除离群值，参数估计的不确定性将降至最低（因为所提出的模型实际上能很好地拟合数据）。
我写了一个非常简单的代码来说明这一点。
当设置 insert_x = [-3] 时，渐近值相对于其量级的标准误差在 print(perr[0]/popt[0]) 附近约为0.36，而在 insert_x = [-30] 附近约为0.014。
如果能提供OP的原始数据，或从共享的图中抓取数据，也能用同样的方法进行演示。
并不需要对原始拟合过程进行专门的微调；所提出的模型相对简单，只包含两个参数；这类情况通常可以通过 curve_fit 来很好地处理。
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
def model(x, asym, k):
return asym-asym*np.exp(-k*x)
x = np.linspace(10, 100, 26)
# define arbitrary number of outliers here
insert_value = [0]
insert_x = [-3]
data = model(x, 300, 0.02)
data = np.insert(data, 0, insert_value)
x = np.insert(x, 0, insert_x)
popt, pcov = curve_fit(model, x, data, p0=[300, 0.1], method='trf')
fit = model(x, *popt)
plt.scatter(x, data, label='Data')
plt.plot(x, fit, c='orange', label='Fitted model ' r"$L_\infty(1-e^{-kt})$")
plt.legend(frameon=False)
plt.show()
perr = np.sqrt(np.diag(pcov))
#print parameter estimates
print(popt)
print(perr)
# SE in the asymptote relative to the magnitude of the asymptote
print(perr[0]/popt[0])
# SE in the saturation constant relative to the magnitude of the saturation constant
print(perr[1]/popt[1])`
看一个绘制的示例：
备选方案2
你正面临两个主要挑战：
- 剪裁过的数据集与当前模型不匹配；
- 离群值影响拟合。
我将认为你的目标是拥有一个能够同时应对这两个挑战、又无需对数据集进行额外处理的流程。
MCVE
首先导入所需的软件包：
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize, stats
我们以剪裁的方式对你的模型进行调整以考虑这一现象：
def model(x, C, k, x0):
return np.clip(C * (1. - np.exp(- k * (x - x0))), 0, None)
我们创建一个合成数据集：
x = np.linspace(-5, 150, 50)
p = (3500, 0.035, 15)
y = model(x, *p)
我们修改一些数据点，以使数据问题与你在图中看到的相同：
y[35] = 0
y[36] = 0
y[40:45] += 50
y[45:50] -= 50
我们首先进行一个朴素的调整，确实会受到离群值的影响。首次调整的目标是捕捉拟合的大致形状，以便在第二步识别离群值并对其进行惩罚。
popt1, pcov1 = optimize.curve_fit(model, x, y)
yhat1 = model(x, *popt1)
我们计算拟合残差的z-score，并构建一个任意权重函数 w 来惩罚离群值：
z = stats.zscore(y - yhat1)
w = (z + 1.) ** 2
我们进行第二次拟合，使用该权重并将 sigma 开关与之一起应用：
popt2, pcov2 = optimize.curve_fit(model, x, y, sigma=w)
# (array([3.48620052e+03, 3.53398594e-02, 1.50469456e+01]),
# array([[ 8.27335558e+02, -2.33783377e-02, -3.84408459e+00],
# [-2.33783377e-02, 1.01213870e-06, 2.51169006e-04],
# [-3.84408459e+00, 2.51169006e-04, 1.11327759e-01]]))
最后对加权模型进行回归：
yhat2 = model(x, *popt2)
图形呈现如下：
正如我们所看到的，第一次拟合受到离群值的影响，而第二次拟合受到的影响要小得多。
相关文章
