2016年4月5日火曜日

[gnuplot] 関数のフィッテイング

gnuplotの機能の一つに,関数のフィッティングがあります.
パラメータを含んだ関数を用意すると,データファイルの数値に対して最適なパラメータの数値を自動的に探索します.
関数は,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 件のコメント :

コメントを投稿