パラメータを含んだ関数を用意すると,データファイルの数値に対して最適なパラメータの数値を自動的に探索します.
関数は,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 件のコメント :
コメントを投稿