12

有限サイズの物理システムのコンピューター シミュレーションを行っています。この後、無限大 (熱力学的限界) への外挿を行っています。一部の理論では、データはシステムのサイズに比例してスケーリングする必要があるため、線形回帰を行っています。

私が持っているデータにはノイズがありますが、データ ポイントごとに誤差範囲を推定できます。たとえば、データ ポイントは次のようになります。

x_list = [0.3333333333333333, 0.2886751345948129, 0.25, 0.23570226039551587, 0.22360679774997896, 0.20412414523193154, 0.2, 0.16666666666666666]
y_list = [0.13250359351851854, 0.12098339583333334, 0.12398501145833334, 0.09152715, 0.11167239583333334, 0.10876248333333333, 0.09814170444444444, 0.08560799305555555]
y_err = [0.003306749165349316, 0.003818446389148108, 0.0056036878203831785, 0.0036635292592592595, 0.0037034897788415424, 0.007576672222222223, 0.002981084130692832, 0.0034913019065973983]

Pythonでこれをやろうとしているとしましょう。

  1. 私が知っている最初の方法は次のとおりです。

    m, c, r_value, p_value, std_err = scipy.stats.linregress(x_list, y_list)

    これにより結果のエラーバーが得られることは理解していますが、これは初期データのエラーバーを考慮していません。

  2. 私が知っている2番目の方法は次のとおりです。 

    m, c = numpy.polynomial.polynomial.polyfit(x_list, y_list, 1, w = [1.0 / ty for ty in y_err], full=False)

ここでは、最小二乗近似で使用される重みとして、各ポイントのエラーバーの逆数を使用します。したがって、ポイントがそれほど信頼できない場合、結果に大きな影響を与えることはありません。これは合理的です。

しかし、これらの両方の方法を組み合わせたものを取得する方法がわかりません。

私が本当に欲しいのは、2番目の方法が何をするかです。つまり、すべてのポイントが異なる重みで結果に影響を与える場合に回帰を使用します。しかし同時に、結果がどれほど正確かを知りたい、つまり、結果の係数の誤差範囲を知りたいのです。

これどうやってするの？

4

4 に答える 4

3

データセットの加重線形回帰を実行する簡潔な関数を作成しました。これは、GSL の「gsl_fit_wlinear」関数を直接翻訳したものです。これは、フィットを実行するときに関数が何をしているかを正確に知りたい場合に役立ちます

def wlinear_fit (x,y,w) :
    """
    Fit (x,y,w) to a linear function, using exact formulae for weighted linear
    regression. This code was translated from the GNU Scientific Library (GSL),
    it is an exact copy of the function gsl_fit_wlinear.
    """
    # compute the weighted means and weighted deviations from the means
    # wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i)
    W = np.sum(w)
    wm_x = np.average(x,weights=w)
    wm_y = np.average(y,weights=w)
    dx = x-wm_x
    dy = y-wm_y
    wm_dx2 = np.average(dx**2,weights=w)
    wm_dxdy = np.average(dx*dy,weights=w)
    # In terms of y = a + b x
    b = wm_dxdy / wm_dx2
    a = wm_y - wm_x*b
    cov_00 = (1.0/W) * (1.0 + wm_x**2/wm_dx2)
    cov_11 = 1.0 / (W*wm_dx2)
    cov_01 = -wm_x / (W*wm_dx2)
    # Compute chi^2 = \sum w_i (y_i - (a + b * x_i))^2
    chi2 = np.sum (w * (y-(a+b*x))**2)
    return a,b,cov_00,cov_11,cov_01,chi2

フィットを実行するには、次のようにします。

a,b,cov_00,cov_11,cov_01,chi2 = wlinear_fit(x_list,y_list,1.0/y_err**2)

aこれは、線形回帰の係数(切片) と(勾配)の最良の推定値とb、共分散行列の要素cov_00およびcov_01を返しcov_11ます。上の誤差の最良の推定値は のa平方根でcov_00あり、上の誤差bは の平方根ですcov_11。残差の加重合計がchi2変数に返されます。

重要: この関数は、データ ポイントの重みとして逆標準偏差ではなく、逆分散を受け入れます。

于 2017-01-11T15:47:04.223 に答える
0

このドキュメントは、独自の加重最小二乗ルーチン (任意のプログラミング言語に適用可能) を理解し、設定するのに役立ちます。

通常、最適化されたルーチンを学習して使用することが最善の方法ですが、ルーチンの内容を理解することが重要な場合もあります。

于 2015-05-25T19:40:48.063 に答える