その他もろもろ (その2)

非線形関数を使った最小自乗フィッティングの方法(その1)．

gnuplot3.7の強力な機能の一つに，関数のフィッティング(当てはめ)があり ます．パラメータを含んだ関数を用意すると，データファイルの数値に対 して最適なパラメータの数値を自動的に探索します．関数はgnuplotで表現 できるものなら，線形・非線形を問いません．ここでは， f(x)=a/(1+b*x*x) というLorentz型の関数を使った例を示します．この関 数のパラメータは，aとbの2 つです．

サンプルのデータとして，a=2, b=1 の値をこの関数に与えて 計算したものを用意しました．この数値に関数をフィットしたときに， a=2, b=1という結果が得られるかどうかを見ようとしてい ます．

-10.000000    0.019802
 -8.947370    0.024674
 -7.894740    0.031582
 -6.842110    0.041828
 -5.789470    0.057941
 -4.736840    0.085333
 -3.684210    0.137236
 -2.631580    0.252359
 -1.578950    0.572561
 -0.526316    1.566160
  0.526316    1.566160
  1.578950    0.572561
  2.631580    0.252359
  3.684210    0.137236
  4.736840    0.085333
  5.789470    0.057941
  6.842110    0.041828
  7.894740    0.031582
  8.947370    0.024674
 10.000000    0.019802

この数値を，"exp.dat" という名前のファイルに保存しておきます．最初 に，gnuplotで当てはめたい関数を定義します．探索するパラメータの変数名 は，xやt等のgnuplotが使う特別な変数でなければ何でも構いません．関数フィッ ティングには， fit コマンドを使います．オプションに，探索した いパラメータ名を via に続けて与えます．

フィットを行うデータファイルに対しては，using, index, every等のオプションが使うことができ，デー タファイル中の任意の数値をフィッティングの対象にすることができます． また，データの誤差を与える事で，各点での重みを変えることができます． データの重みは，誤差σの2乗の逆数(σ^{-2})となります．重み付き最小 自乗法を行うためには，データファイルの3カラム目に誤差を用意し， fit f(x) "exp.dat" using 1:2:3 via a,b のようにusing を使って誤差のカラムを指定します．誤差を与えない場合は，全データの 重みは1になります．

gnuplot> f(x)=a/(1+b*x*x)  
gnuplot> fit f(x) "exp.dat" via a,b


Iteration 0
WSSR        : 1.43879           delta(WSSR)/WSSR   : 0
delta(WSSR) : 0                 limit for stopping : 1e-05
lambda	  : 0.201184

                         .......

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a               = 2                +/- 5.236e-07    (2.618e-05%)
b               = 0.999998         +/- 7.397e-07    (7.397e-05%)


correlation matrix of the fit parameters:

               a      b      
a               1.000 
b               0.837  1.000 

途中，フィッティングのログが長々と画面を流れて行きます．同じ内容が， fit.log という名前のファイルにも同時に書き出されています．フィッティ ングがうまく終れば，得られたパラメータの値と，それらの相関(共分散) 行列が出力されます．上の例では，a=2, b=1がほぼ正しく 得られています．

得られた数値は，変数a,bに入っていますので，関数f(x)をデータと ともにプロットすると，フィッティングの結果をプロットできます．

gnuplot> plot f(x) with lines, "exp.dat" using 1:2 w points

非線形関数を使った最小自乗フィッティングの方法(その2)．

次はもう少し具体的な例を用います．次のような実験データがあります．こ のデータに関数f(x)=c*(x-a)**bをフィットし，変数a,b,cの値を求 めます．

# X        Y     Z
  6.295    8.4   0.3
  6.795   23.9   0.7
  7.826  104.0   3.0
  8.830  255.0   8.0
  9.841  421.0  13.0
 10.860  566.0  18.0
 11.864  690.0  22.0

各データ点の重みは，誤差の値をZとすると，1/Z*Zとなります．従って誤 差が大きい程重みは小さくなります．

まず，実験値だけをプロットしてその傾向を見ます．

gnuplot> set log y
gnuplot> plot "exp.dat" using 1:2:3 w yerrorbars

次に変数a,b,cの初期値を決めます．gnuplotのフィッティングでは， 初期値を与えていないと全て1として計算を始めます．しかし，反復回数が 多くなったり場合によっては結果が出なかったりしますので，ほどほどの 初期値を決めておいた方がうまく行きます．f(x)=c*(x-a)**bをあてはめま すので，上のプロットを見れば，aは大体6位だろうと想像できます (x=aで関数はゼロになるはずなので)．xが1増える と数十倍になっているので，bは1〜2程度でしょう．また， cは10のオーダー(f(x)の対数をとったら分かります)，これを初期 値にして，fit を実行します．

gnuplot> a=6
gnuplot> b=1
gnuplot> c=10
gnuplot> f(x)=c*(x-a)**b
gnuplot> fit f(x) "exp.dat" using 1:2:3 via a,b,c

17回の反復の後に計算が収束したメッセージが出，最終的なパラメータの値 とその共分散が出力されます．得られた値はa=5.77, b=1.89, c=25.8 でしたので，仮定した初期値はまずまずだっ たようです(プロットしてみるとあまりデータにあっていないのですが．．．)

After 17 iterations the fit converged.
final sum of squares of residuals : 90.8196
rel. change during last iteration : -9.71916e-09

degrees of freedom (ndf) : 4
rms of residuals      (stdfit) = sqrt(WSSR/ndf)      : 4.76497
variance of residuals (reduced chisquare) = WSSR/ndf : 22.7049

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a               = 5.76731          +/- 0.176        (3.051%)
b               = 1.89369          +/- 0.2127       (11.23%)
c               = 25.7956          +/- 9.807        (38.02%)


correlation matrix of the fit parameters:

               a      b      c      
a               1.000 
b              -0.944  1.000 
c               0.975 -0.975  1.000 

この結果を，先程のプロットに重ねて表示します．

gnuplot> plot f(x),"exp.dat" using 1:2:3 w yerrorbars