在 lmfit 中考虑实验误差

Question

我正在尝试将lmfit 实施到我的拟合例程中，但在定义错误时遇到了问题。我假设我在这个平台上阅读了有关该主题的先前问题，并且我也浏览了文档，但我的一些疑问仍然存在。

下面是我想要实现的一个完整且最小的示例。

import corner
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import scipy as sp
import lmfit


font = {'fontname':'candara', "fontweight":"light"}
plt.style.use('ggplot')
plt.rcParams["figure.figsize"] = (8,4)
ax_fit_kws = dict(xlim=(0,0.12), ylim=(0.5,1.1))
ax_res_kws = dict(xlim=(0,0.12), ylim=(-0.1,0.1))

def mono_exp(SL_array, a, b):
    model = a * np.exp(-b * SL_array)
    return model
model = lmfit.Model(mono_exp)

SL_array = np.array((0.030, 0.040, 0.060, 0.080, 0.10))
data= np.array((1., 0.9524336, 0.92452666, 0.87995659, 0.82845576))
errs = np.array((0.00029904, 0.00049384, 0.00076344, 0.00053886, 0.00066012))

params = model.make_params(a=0, b=0)

result = model.fit(data=data, params=params, SL_array=SL_array, method="Nelder", markersize=10, weights=errs)

lmfit.report_fit(result)
result.plot(yerr = errs, ax_fit_kws=ax_fit_kws, ax_res_kws=ax_res_kws)

emcee_kws = dict(steps=400, burn=30, thin=20, is_weighted=False,
                 progress=True)
emcee_params = result.params.copy()
emcee_params.add('__lnsigma', value=np.log(0.1), min=np.log(0.001), max=np.log(2.0))
result_emcee = model.fit(data=data, SL_array=SL_array, params=emcee_params, method='emcee',
                         nan_policy='omit', fit_kws=emcee_kws)

lmfit.report_fit(result_emcee)

ax = plt.plot(SL_array, model.eval(params=result.params, SL_array=SL_array), label='Nelder', zorder=100)
result_emcee.plot_fit(ax=ax, data_kws=dict(color='gray', markersize=10), yerr = errs)
emcee_corner = corner.corner(result_emcee.flatchain, labels=result_emcee.var_names,
                             truths=list(result_emcee.params.valuesdict().values()))
plt.show()

我的问题相当简单：我希望最初的Nelder 拟合例程将errs 数组作为data 数组上的错误考虑（这是我通过实验确定的点）。我不确定调用weights=errs 能否实现这一目标。我已经尝试过这里实现的解决方案：How do I include errors for my data in the lmfit least squares miniimization, and what is this error for conf_interval2d function in lmfit? 但我无法让它工作。

我不太清楚的另一点是：我适合的emcee 部分是否考虑了Nelder 例程的残差？

非常感谢！

编辑

经过更多研究，我现在认为，在调用权重时，我实际上应该给出1/err。通过实施此更改，并应用 scale_covar=False 如别处所述（How to properly get the errors in lmfit），同时人为地增加错误值（例如，故意将 errs 数组乘以 100 倍）我确实得到了拟合参数错误增加实质上，这是一种预期的行为。简而言之：result = model.fit(data=data, params=params, SL_array=SL_array, method="Nelder", markersize=10, weights=errs) 已更改为 result = model.fit(data=data, params=params, SL_array=SL_array, method="Nelder", markersize=10, weights=1/errs)。这是正确的吗？

在我的情况下，我仍然对 emcee 的实现感到困惑。

Answer 1

您在编辑中所说的内容是正确的：您想使用 weights=1./err 通过数据中的不确定性来正确加权 data 和 model 的残差 err。

您可能也想在对model.fit(..., method='emcee') 的调用中使用相同的内容。

我应该说emcee 在lmfit 中的使用相当令人困惑，并且给人的不幸印象是它很合适。这根本不是真的，因为emcee（以及真正的 MCMC 作为一种方法）在“系统地优化参数值以找到改进的解决方案”的意义上并不能真正做到合适。它所做的是探索输入参数值附近的参数空间（这恰好是Nelder 方法的解决方案）。 这种探索可能会找到（更像是“偶然发现”而不是“寻找”）改进的解决方案，并且结果将反映它所做的探索。