gnuplot

2017年3月11日 (土)

[gnuplot]重み付き最小二乗法による回帰モデルを作成する(2変量、目的変数にのみ誤差あり)

いわゆる y = f(x) のような方程式で、従属変数(この場合は変数y)にのみ誤差がある場合の、重み付き最小二乗法をfitコマンドで行う方法。

以下のようにデータファイルを作成して(ここでファイル名はxye.txtとする)、fitコマンドのusingオプションを以下のように使用すればよい。なお、fitコマンドにおける重みの指定は、重みではなく誤差を指定することに注意(Rのlm関数とは違う)。fitコマンドは指定された誤差を対応する変数の標準偏差σと見なし、各データの重みを1/σ^2として計算をする。

以下、計算例のためのデータ。左から順に、独立変数(x)の値、従属変数(y)の値、誤差(e)の値。

1  1  0.5
2  4  1.0
3  9  0.5
5 25 10.0

実行例。

gnuplot> f(x) = a + b * x                              
gnuplot> fit f(x) "xye.txt" using 1:2:3 via a,b
         warning:
        > Implied independent variable y not found in fit function.
        > Assuming version 4 syntax with zerror in column 3 but no zerror keyword.

iter      chisq       delta/lim  lambda   a             b            

   0 1.5387570311e+00  0.00e+00 1.01e+01 -3.162299e+00  4.030045e+00
   * 1.5387570311e+00  1.15e-10 1.01e+02 -3.162299e+00  4.030045e+00
   * 1.5387570311e+00  4.33e-11 1.01e+03 -3.162299e+00  4.030045e+00
   * 1.5387570311e+00  4.33e-11 1.01e+04 -3.162299e+00  4.030045e+00
   1 1.5387570311e+00  0.00e+00 1.01e+03 -3.162299e+00  4.030045e+00
iter      chisq       delta/lim  lambda   a             b            

After 1 iterations the fit converged.

final sum of squares of residuals : 1.53876
rel. change during last iteration : 0

degrees of freedom    (FIT_NDF)                        : 2

rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.877142
variance of residuals (reduced chisquare) = WSSR/ndf   : 0.769379
p-value of the Chisq distribution (FIT_P)              : 0.463301

Final set of parameters            Asymptotic Standard Error

=======================            ==========================
a               = -3.1623          +/- 0.6834       (21.61%)
b               = 4.03005          +/- 0.3084       (7.652%)

correlation matrix of the fit parameters:

                a      b      
a               1.000
b              -0.904  1.000

図を表示する。

gnuplot> plot [0:6][0:30] "xye.txt" with points pointtype 7 
gnuplot> replot f(x)

Figure11

紫丸が元データ。緑直線が重みを考慮したときの回帰直線。データを見ての通り、4番目の組(x = 5)は誤差を大きくしているので、重みを考慮させると、残り3点で決められたような回帰モデルになっているのがわかる。