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

Question

我刚开始学习Python，想用lmfit这个包来检查我自己的计算。不过我有两个问题：第一，我不知道怎么把数据的误差（sig）加进去；第二，我在使用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）加进去，这样我才能计算加权和简化的卡方值呢？

第二部分：此外，当我尝试使用以下代码来绘制置信区间时，

xs, ys, grid = conf_interval2d(mi, 'a', 'b', 20, 20)

我遇到了以下错误：

* ValueError: 创建意图（cache|hide)|optional数组失败——必须定义维度，但得到了(0,)

或者

例程DGESV的参数4不正确

Mac OS BLAS参数错误在DGESV中，参数#0，（不可用），是0

Answer 1

在residual()函数中，你需要自己根据误差来给数据加权。

根据lmfit的文档（虽然不太容易找到）：

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

不过，这并不难。根据维基百科关于加权最小二乘法拟合的介绍：

如果测量值之间没有相关性，但不确定性不同，可以采用一种修改的方法。Aitken证明，当最小化加权的平方残差和时，如果每个权重等于测量方差的倒数，那么结果是最佳线性无偏估计（BLUE）。

但是，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

这样应该能给你正确加权的拟合结果。