4

我是 python 新手，并试图使用 lmfit 包来检查我自己的计算，但是我不确定（1）如何为以下测试（和 2）包含数据错误（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）的错误，以便我可以计算加权和减少的卡方？

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

我收到以下错误：

* ValueError: failed to create intent(cache|hide)|optional array-- 必须有定义的维度但是得到 (0,)

或者

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

4

1 回答 1

5

您必须自己根据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

这应该为您提供正确加权的拟合。

于 2013-04-12T19:35:01.823 回答