在 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 作为一种方法）不能真正做到 "systematically refine parameter values in order to find an improved solution" 意义上的拟合。它正在做的是探索输入参数值附近的参数 space（恰好是 Nelder 方法的解）。此探索可能会发现（更像是 "stumble upon" 而不是 "seek"）改进的解决方案，结果将反映它所做的探索。

在 lmfit 中考虑实验误差

Taking experimental errors into account in lmfit

curve-fitting

python-3.x

lmfit