前回はデータの補間線について紹介したが、今回はデータに対して指定の関数にフィッティングする方法を紹介する。実験データなどにある既知の工学的なモデル（関数）をフィッティングすることで、実験データがモデルに従うか、また得られる係数からどのような特徴があるかなど評価する上でデータのフィッティングはよく行われている。 atatat.hatenablog.com
1. 線形フィッティング
フィッティングする関数が線形であれば線形フィッティング、非線形であれば非線形フィッティングと呼ばれているが、まず線形フィッティングについて紹介する。データファイルlinea_dat.txtはExcelでXは0～15,Yは2*X+3+(RAND()-0.5)*5としてデータを生成した。勾配2, 切片3の線形データに乱数として-2.5～2.5を発生させてノイズとして加えている。
#X Y linear_dat.txt 0 3.701591018 1 6.724659473 2 4.98067816 3 8.49673559 4 10.90453657 5 14.01262334 6 14.31811822 7 17.74440048 8 18.05810669 9 20.55938899 10 25.49320738 11 22.67537007 12 27.61579479 13 26.82399374 14 31.01257974 15 33.18377999
このデータをプロットすると以下のようになる。
gnuplotを用いて関数にフィッティングを行うためには、フィッティングする関数とパラメータを定義する。今回は線形関数なのでax+bとしてフィッティングにより係数a,bを求める。 フィッティング関数f(x)を定義するには
f(x)=a*x+b
とする。次にデータファイルlinear_dat.txtをf(x)にフィッティングして、係数a,bを算出するにはfit 定義した関数 データファイル名 via パラメータ1,パラメータ2,...とする。
fit f(x) 'linear_dat.txt' via a,b
とするとコンソール上に以下のように出力される。
iter chisq delta/lim lambda a b
0 2.7830299416e+01 0.00e+00 1.27e+01 2.008065e+00 2.882490e+00
1 2.7193134224e+01 -2.34e+03 1.27e+00 1.983556e+00 2.943640e+00
2 2.6822558511e+01 -1.38e+03 1.27e-01 1.955864e+00 3.221762e+00
3 2.6821874903e+01 -2.55e+00 1.27e-02 1.954648e+00 3.234238e+00
4 2.6821874903e+01 -5.13e-07 1.27e-03 1.954647e+00 3.234244e+00
iter chisq delta/lim lambda a b
After 4 iterations the fit converged.
final sum of squares of residuals : 26.8219
rel. change during last iteration : -5.12538e-12
degrees of freedom (FIT_NDF) : 14
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 1.38414
variance of residuals (reduced chisquare) = WSSR/ndf : 1.91585
Final set of parameters Asymptotic Standard Error
======================= ==========================
a = 1.95465 +/- 0.07507 (3.84%)
b = 3.23424 +/- 0.6608 (20.43%)
correlation matrix of the fit parameters:
a b
a 1.000
b -0.852 1.000
フィッティングに関する様々な情報が表示されるが、基本的にはフィッティングが収束したかどうか:After 4 iterations the fit converged.および算出した係数の値:Final set of parametersが重要な情報となる。
今回は4回の繰り返し計算で収束し、係数a=1.95465 , b=3.23424というのが結果となる。ファイルデータ生成時にa=2, b=3としているのでノイズの影響はあるがある程度正しい値が算出できている。フィッティングした関数f(x)がデータにうまく近似できているかどうか確認するために、linear_dat.txtとf(x)をプロットする。
plot 'linear_dat.txt', f(x)
係数a,bについてはフィッティング時にa=1.95465 , b=3.23424が割り当てられているの、そのままf(x)をプロットすればよい。 以下がプロットした結果で、よく近似できていることが確認できる。
2. 非線形フィッティング
次に非線形フィッティングについて実行する。データファイルexp_dat.txtはExcelでXは0～15,Yは1.5*EXP(0.25*X)*(1+0.3*(RAND()-0.5))としてデータを生成した。指数関数に対してその値の3割がノイズになるように生成している。（データのプロットは省略）
#X Y exp_dat.txt 0 1.571329315 1 1.893034492 2 2.495464388 3 3.582999282 4 3.875627378 5 5.689683543 6 7.462634119 7 9.534954001 8 10.79657878 9 14.78070029 10 18.37349399 11 26.01796549 12 32.61328597 13 43.84045216 14 46.07542661 15 55.27057992
今回はフィッティング関数をc*exp(d*x)としてフィッティングにより係数c,dを求める。非線形関数であっても線形関数と同様に行う。
g(x)=c*exp(d*x) fit g(x) 'exp_dat.txt' via c,d
以下実行すると、フィッティングの途中結果が大量に出力されて、最後に以下のようになる。
After 102 iterations the fit converged.
final sum of squares of residuals : 67.7492
rel. change during last iteration : -1.29457e-06
degrees of freedom (FIT_NDF) : 14
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 2.19983
variance of residuals (reduced chisquare) = WSSR/ndf : 4.83923
Final set of parameters Asymptotic Standard Error
======================= ==========================
c = 2.21331 +/- 0.2982 (13.47%)
d = 0.218301 +/- 0.01006 (4.607%)
correlation matrix of the fit parameters:
c d
c 1.000
d -0.986 1.000
繰り返し計算を102回行って、係数c=2.21331 , d= 0.218301 という結果が得られた。ファイルデータ生成時はc=1.5, b=0.25であったのでそこそこ正解に近いといえる。
実際にフィッティングした関数g(x)とexp_dat.txtをプロットしてみると、よく近似できていることがわかる。
3. 線形フィッティングとの違い
手順自体は線形フィッティングと同じだが繰り返し回数が明らかに多く、非線形関数である指数関数の方が収束しにくい。そこである程度係数c,dの値がわかっている場合は近い値を初期値として与えることで、早く収束することがある。例えばc=3,d=0.1とすると
c=3 d=0.1 fit g(x) 'exp_dat.txt' via c,d
実行すると以下のように繰り返し回数が6回に減り、c,dもほぼ同じ値に収束する。
iter chisq delta/lim lambda c d 0 4.9785030725e+03 0.00e+00 8.44e+00 3.000000e+00 1.000000e-01 * 4.2978399981e+04 8.84e+04 8.44e+01 2.690339e+00 2.876011e-01 1 2.6596547584e+03 -8.72e+04 8.44e+00 3.596992e+00 1.271493e-01 2 7.0017603817e+02 -2.80e+05 8.44e-01 2.239844e+00 2.343823e-01 3 7.3707863613e+01 -8.50e+05 8.44e-02 2.220472e+00 2.198945e-01 4 6.7751037456e+01 -8.79e+03 8.44e-03 2.217328e+00 2.181890e-01 5 6.7749239406e+01 -2.65e+00 8.44e-04 2.213105e+00 2.183081e-01 * 6.7749244663e+01 7.76e-03 8.44e-03 2.213427e+00 2.182973e-01 * 6.7749244663e+01 7.76e-03 8.44e-02 2.213427e+00 2.182973e-01 * 6.7749244662e+01 7.76e-03 8.44e-01 2.213427e+00 2.182973e-01 * 6.7749244629e+01 7.71e-03 8.44e+00 2.213427e+00 2.182973e-01 * 6.7749242373e+01 4.38e-03 8.44e+01 2.213355e+00 2.182998e-01 6 6.7749239285e+01 -1.78e-04 8.44e+00 2.213117e+00 2.183078e-01 iter chisq delta/lim lambda c d After 6 iterations the fit converged. final sum of squares of residuals : 67.7492 rel. change during last iteration : -1.78413e-09 degrees of freedom (FIT_NDF) : 14 rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 2.19983 variance of residuals (reduced chisquare) = WSSR/ndf : 4.83923 Final set of parameters Asymptotic Standard Error ======================= ========================== c = 2.21312 +/- 0.2981 (13.47%) d = 0.218308 +/- 0.01006 (4.607%)
ただし、はずれた値を初期値として指定すると、例えばc=10,d=9とすると
iter chisq delta/lim lambda c d 0 1.8176494128e+119 0.00e+00 1.09e+61 1.000000e+01 9.000000e+00 1 2.9703912914e+118 -5.12e+05 1.09e+60 9.999536e+00 8.939622e+00 iter chisq delta/lim lambda c d The maximum lambda = 1.000000e+20 was exceeded. Fit stopped. final sum of squares of residuals : 2.97039e+118 rel. change during last iteration : -5.11923 degrees of freedom (FIT_NDF) : 14 rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 4.6062e+58 variance of residuals (reduced chisquare) = WSSR/ndf : 2.12171e+117 Final set of parameters Asymptotic Standard Error ======================= ========================== c = 9.99954 +/- 2.875e+05 (2.875e+06%) d = 8.93962 +/- 1791 (2.003e+04%)
"The maximum lambda = 1.000000e+20 was exceeded. Fit stopped."と表示されてフィッティングが打ち切られたことがわかる。 非線形関数でフィッティングする場合は、ある程度近い初期値を指定することが重要になる。
4. まとめ
今回はgnuplotを用いた線形フィッティングおよび非線形フィッティングについて紹介した。非線形フィッティングに関しては、線形フィッティングと同様の手順で実行できるが適切な初期値を与えないと、うまく収束しないことを取り上げた。
