我如何在 lmfit 最小二乘最小化中包含我的数据的错误，以及 lmfit 中 conf_interval2d 函数的错误是什么？

Question

我是 python 的新手，并尝试使用 lmfit 包来检查我自己的计算，但是我不确定 (1) 如何包含数据 (sig) 的错误以用于以下测试（和 2）的错误我使用如下所示的 conf_interval2d 获取）：

    import numpy as np
    from lmfit import Parameters, Minimizer, conf_interval, conf_interval2d, minimize, printfuncs


    x=np.array([ 0.18,  0.26,  1.14,  0.63,  0.3 ,  0.22,  1.16,  0.62,  0.84,0.44,  1.24,  0.89,  1.2 ,  0.62,  0.86,  0.45,  1.17,  0.59, 0.85,  0.44])
    data=np.array([ 68.59,  71.83,  22.52,44.587,67.474 ,  55.765,  20.9,41.33783784,45.79 ,  47.88,   6.935,  34.15957447,44.175,  45.89230769,  57.29230769,  60.8,24.24335594,  34.09121287,  42.21504003,  26.61161674])
    sig=np.array([ 11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409])

    def residual(pars, x, data=None):
        a=pars['a'].value
        b=pars['b'].value
        model = a + (b*x)
        if data is None:
            return model
        return model-data

    params=Parameters()
    params.add('a', value=70.0)
    params.add('b', value=40.0)

    mi=minimize(residual, params, args=(x, data))
    #mi=minimize(residual, params, args=(x,), kws={'data': data})#is this more correct?
    ci, trace = conf_interval(mi, trace=True)

到目前为止，这工作正常，但如上所述，我如何包括数据 (sig_chla) 的错误，以便我可以计算加权和减少的卡方？

第 2 部分：进一步，当我尝试使用以下内容以便我可以绘制置信区间 xs, ys, grid = conf_interval2d(mi, 'a', 'b', 20, 20)

我收到以下错误：

* ValueError：未能创建意图（缓存|隐藏）|可选数组--必须定义维度但得到（0，）

或者

例程 DGESV 的参数 4 不正确 DGESV 中的 Mac OS BLAS 参数错误，参数 #0（不可用）为 0

Answer 1

您必须自己根据residual()函数中的错误对数据进行加权。

从lmfit文档（虽然不是很容易找到）：

请注意，卡方和缩减卡方的计算假设返回的残差函数已根据数据中的不确定性进行了适当的缩放。为了使这些统计数据有意义，编写要最小化的函数的人必须适当地缩放它们。

不过，这并不难。来自维基百科关于加权最小二乘拟合的条目：

但是，如果测量结果不相关但具有不同的不确定性，则可能会采用修改后的方法。Aitken 表明，当残差的加权平方和最小化时，如果每个权重都等于测量方差的倒数，则为蓝色。

然而，lmfit接受残差，而不是平方残差，所以而不是仅仅去

    # This is what you do with no errorbars -> equal weights.
    resids = model - data
    return resids

你应该做这样的事情（scipy导入为sp）：

    # Do this to include errors as weights.
    resids = model - data
    weighted = sp.sqrt(resids ** 2 / sig ** 2)
    return weighted

这应该会给你正确的加权合身。