パラメータを含んだ関数を用意すると，データファイルの数値に対して最適なパラメータの数値を自動的に探索します．
関数は，gnuplotで表現できるものなら，線形・非線形を問いません．
具体例を示します．
以下のようなデータ（圧密試験結果，data.csv）があったとします．
# Consolidation Pressure p(kN/m2),Void Ratio e
0,0.593
9.8,0.584
19.6,0.58
39.2,0.573
78.5,0.564
157,0.551
314,0.536
628,0.515
1256,0.484
このデータに関して，関数
f(x)=c*(x-a)**b
をフィッティングして，変数a, b, cの値を求めます．
まずは，以下のようなスクリプトを実行してデータをプロットしてみます．
set logscale x 10
set xrange [1:10000]
set tics font "Times New Roman,25"
set yrange [0.35:0.65]
set ytics 0.05
set mytics 5
set xlabel "{/=30 Consolidation stress p(kN/m^{2})}"
set ylabel "{/=30 Void Ratio e}"
set size square
set key font "Times New Roman,25"
set key right top
plot "./data.csv" using 1:2 w p ps 1.5 title "Sample 1"
次に変数a,b,cの初期値を求めます．
gnuplotのフィッティングでは，初期値を与えないと全て"1"として計算を始めます．
しかし，反復回数が多くなったり，結果が得られないこともあるようなので，当たりをつけて初期値を決めた方が良いようです（以下の例では初期値を与えない場合を示します）．
最初に，フィッテイングする関数を定義します．
探索するパラメータの変数名は，xやt等のgnuplotが使う特別な変数でなければ何でも大丈夫です．
関数フィッティングには，fitコマンドを使います．
オプションに，探索したいパラメータ名を via に続けて与えます．
f(x)=c*(x-a)**b
fit f(x) "./data.csv" using 1:2 via a,b,c
フィットを行うデータファイルに対しては，using, index, every等のオプションが使うことができ，データファイル中の任意の数値をフィッティングの対象にできます．
また，データの誤差を与える事で，各点での重みを変えることができます．
データの重みは，誤差σの2乗の逆数(σ^{-2})となります．
重み付き最小 自乗法を行うためには，データファイルの3カラム目に誤差を用意し， fit f(x) "exp.dat" using 1:2:3 via a,b のようにusing を使って誤差のカラムを指定します．誤差を与えない場合は，全データの 重みは1になります．
以下のスクリプト（test.gp）を実行すると，下図が得られます．
set datafile separator ","
set logscale x 10
set xrange [1:10000]
set tics font "Times New Roman,25"
set yrange [0.35:0.65]
set ytics 0.05
set mytics 5
set xlabel "{/=30 Consolidation stress p(kN/m^{2})}"
set ylabel "{/=30 Void Ratio e}"
set size square
set key font "Times New Roman,25"
set key right top
plot "./Dc90-2.csv" using 1:2 w p ps 1.5 title "Sample 1"
f(x)=c*(x-a)**b
fit f(x) "./Dc90-2.csv" using 1:2 via a,b,c
replot f(x), "./Dc90-2.csv" using 1:2 w yerrorbars
ちなみに，Mac OS Xのターミナルから，上記のスクリプトを実行した場合，ターミナルには以下のように表示され，イタレーションの回数（今回の例では15回）や，フィッテングされた変数の値（a, b, c )などがわかります．
gnuplot> load 'test.gp'
iter chisq delta/lim lambda a b c
0 2.0958191547e+06 0.00e+00 2.07e+03 1.000000e+00 1.000000e+00 1.000000e+00
1 3.1401236259e+05 -5.67e+05 2.07e+02 1.000059e+00 8.652118e-01 9.785288e-01
2 4.3661665492e+04 -6.19e+05 2.07e+01 1.000868e+00 7.302030e-01 9.158026e-01
3 3.9742713760e+03 -9.99e+05 2.07e+00 1.008882e+00 6.584271e-01 4.525282e-01
4 7.1730628717e+01 -5.44e+06 2.07e-01 9.923941e-01 6.351242e-01 7.812126e-02
5 8.6890175128e+00 -7.26e+05 2.07e-02 6.980505e-01 4.985268e-01 7.998705e-02
6 5.2167471432e-01 -1.57e+06 2.07e-03 -7.175606e+00 2.711390e-01 1.407226e-01
7 3.0807202872e-01 -6.93e+04 2.07e-04 -2.955954e+01 2.491715e-02 3.138176e-01
8 1.1567525775e-01 -1.66e+05 2.07e-05 -1.338551e+01 -9.192671e-02 6.734814e-01
9 5.8730090513e-03 -1.87e+06 2.07e-06 -3.335614e+01 -4.813472e-02 7.368142e-01
10 9.9947209776e-05 -5.78e+06 2.07e-07 -8.772884e+01 -6.842443e-02 7.982868e-01
11 8.0106216817e-05 -2.48e+04 2.07e-08 -1.408457e+02 -8.119137e-02 8.727729e-01
12 3.6356267939e-05 -1.20e+05 2.07e-09 -1.692429e+02 -8.667967e-02 9.118675e-01
13 3.3880317042e-05 -7.31e+03 2.07e-10 -1.775257e+02 -8.825277e-02 9.234001e-01
14 3.3854660121e-05 -7.58e+01 2.07e-11 -1.791170e+02 -8.857467e-02 9.256578e-01
15 3.3854399746e-05 -7.69e-01 2.07e-12 -1.793814e+02 -8.862971e-02 9.260331e-01
iter chisq delta/lim lambda a b c
After 15 iterations the fit converged.
final sum of squares of residuals : 3.38544e-05
rel. change during last iteration : -7.69104e-06
degrees of freedom (FIT_NDF) : 5
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 0.00260209
variance of residuals (reduced chisquare) = WSSR/ndf : 6.77088e-06
Final set of parameters Asymptotic Standard Error
======================= ==========================
a = -179.381 +/- 47.19 (26.3%)
b = -0.0886297 +/- 0.01025 (11.57%)
c = 0.926033 +/- 0.06838 (7.384%)
correlation matrix of the fit parameters:
a b c
a 1.000
b 0.968 1.000
c -0.979 -0.999 1.000
gnuplot>
0 件のコメント :
コメントを投稿