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 乗を計算できるようにするにはどうすればよいでしょうか?

パート 2: さらに、信頼区間 xs, ys, grid = conf_interval2d(mi, 'a', 'b', 20, 20) をプロットできるように、次を使用しようとすると

次のエラーが表示されます。

* ValueError: インテント (キャッシュ|非表示) | オプションの配列の作成に失敗しました -- 次元を定義する必要がありますが、(0,) を取得しました

また

ルーチン DGESV のパラメータ 4 が正しくありませんでした DGESV の Mac OS BLAS パラメータ エラー、パラメータ #0、(使用不可)、0 です

4

1 に答える 1

5

関数内で自分でエラーによってデータを重み付けする必要がありますresidual()。

ドキュメントからlmfit（ただし、見つけるのは簡単ではありません）：

カイ 2 乗と簡約カイ 2 乗の計算では、返された残差関数がデータの不確実性に対して適切にスケーリングされていることを前提としていることに注意してください。これらの統計を意味のあるものにするには、最小化する関数を作成する人が統計を適切にスケーリングする必要があります。

とはいえ、それほど難しいことではありません。加重最小二乗フィッティングに関するウィキペディアのエントリから:

ただし、測定値に相関がないが不確実性が異なる場合は、修正されたアプローチが採用される可能性があります。Aitken は、残差の二乗の加重和が最小化されたときに、各加重が測定値の分散の逆数に等しい場合、青色であることを示しました。

ただし、lmfit残差の二乗ではなく、残差を取ります。

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

次のようにする必要があります（scipyimport as を使用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 に答える