gnuplot / intro / plotfunc

gnuplot 入門 — 関数表示編

関数の定義

式で書けるような関数であれば,gnuplotはその式を計算して表示すること
ができます.「式で書ける」というのは,例えばf(x)=a*x+b のような形で
書けることを表し,微分や積分などを含むものはgnuplotだけでは計算でき
ません.gnuplotは三角関数やBessel関数,Gamma関数等,多くの関数を内
部に持っており,それら自身,あるいはそれらを組み合わせたグラフを表
示することができます.

関数のグラフそのものを利用することは(数学の授業目的以外には)少ないと
思います(何かあるかもしれません.あったら教えてください).gnuplotの
関数機能の威力が最も発揮されるのは,データへの関数当てはめ(データ
フィッティング),つまり最小自乗法です.実験データを,直線等の単純な
解析曲線で近似することがしばしばありますが,gnuplotを使えば非線形関
数の当てはめも簡単に行うことができます.この章では関数の定義の方法
と,その関数に含まれるパラメータに対する最小自乗法を説明します.

次式で表される,Lorentz型関数と 1/sqrt(x) の和の関数を表示してみます.
式に含まれるa,b,c,dの4つの変数は,実験データから求められるものとし
ます.

gnuplot / intro / plotfunc

gnuplotで関数を定義するには,FortranやCなどのプログラミングと全く同
じように式を書きます.下の例の
f(x)=c/((x-a)*(x-a)+b)+d/sqrt(x) が,上の関数を定義してい
る部分です. (x-a)の2自乗の部分は,Fortran風に
(x-a)**2
と書くこともできます.変数 a,b,c,d は任意ですが,とり
あえず適当な数値を与えています.

 gnuplot> a=0.25 gnuplot> b=0.02 gnuplot> c=0.05 gnuplot> d=0.1 gnuplot> f(x)=c/((x-a)*(x-a)+b)+d/sqrt(x) gnuplot> set xrange [0:1] gnuplot> set yrange [0:4] gnuplot> plot f(x) 
plotfunc1
up

関数値を求める

gnuplotはユーザが定義した式を内部で計算しますが,グラフを表示するだけでなく,
数値を直接表示することもできます.あるX座標の値を求めるには,print
コマンドを用います.式に含まれるパラメータ(a,b,c,d)を変えれば,関数値も
変化します.なお,計算は倍精度で行われます.

 gnuplot> print f(0.25) 2.7 gnuplot> print f(0.4) 1.33458447124371 gnuplot> a=0.4 gnuplot> print f(0.4) 2.65811388300842 

計算された関数値をスプレッドシート等に取り込んで利用したいこともあり
ます.数値を表にして出力したいなら,特殊なterminal, table
を用います. set output でファイルを指定すれば,
計算結果はそのファイルに書き込まれます.

 gnuplot> set term table Terminal type set to 'table' gnuplot> plot f(x) #Curve 0, 100 points #x y type 0 0 u 0.010101 1.63972 i 0.020202 1.39031 i 0.030303 1.30688 i ....  0.979798 0.191506 i 0.989899 0.188622 i 1 0.185837 i  gnuplot> set output "calc.plt" gnuplot> replot 
up

最小自乗法でパラメータを求める

パラメータ a,b,c,dに対し,実験データに基づいて最適値を求めます.データ
はファイル exp.dat に用意しておきます.

 2.5000E-03 3.0420E+00 6.47E-01 3.5000E-03 2.5700E+00 4.37E-01 4.5000E-03 2.3020E+00 2.53E-01 ...  7.0000E-01 2.7420E-01 2.14E-03 7.5000E-01 2.5680E-01 1.81E-03 8.0000E-01 2.4630E-01 1.59E-03 

データファイルには,(x,y,z)の組が1行に書かれています.X,Y座標に加え
て,Yの誤差Zが3カラム目に与えられています.誤差は絶対誤差で,データ
Yと同じ単位を持ちます.例えば,上の1行目で,Yが3.04cmなら,誤差は
0.647cmといった具合です.誤差の逆数が,各データ点の重みとなります.
誤差が無ければ,全てのデータは同じ重みを持つものと解釈されます.

まずはデータと関数を同時に重ねてプロットします.
実験データ編
数値計算編で述べたように,
軸名等を適宜調整しています.

 gnuplot> set xlabel "Energy [MeV]" gnuplot> set ylabel "Cross Section [b]" gnuplot> set xtics 0.1 gnuplot> set ytics 0.5 gnuplot> plot f(x) title "Lorentzan", > "exp.dat" using 1:2:3 title "experiment" with yerrors 

plotfunc2

上で,パラメータ a,b,c,d は適当に与えたと書きましたが,実はこのデー
タに対して大体の数値を決めたものです.aはLorentz型のピーク位置に相
当しますが,この実験データの図から,大体0.25としました.bの平方根は
ピークの幅に対応します.図から幅は0.1程度,その平方が0.01ですが,
ちょっと大きめにして0.02を入れてみました.


gnuplotを使って最小自乗法をするのはとても簡単です. fit
マンドを用い,求めたいパラメータを via というオプションに
続けて書くだけです.但し,この関数のような非線形関数を使って多くの
パラメータを同時にフィットには,初期値に注意する必要があります.初
期値が最適値から非常に離れている場合は,最小自乗解を求めるのが難し
くなります.上の数値をそのまま使って最小自乗法を適用すると,次のよ
うになります.

 gnuplot> fit f(x) "exp.dat" using 1:2:3 via a,b,c,d   Iteration 0 WSSR        : 96618.1           delta(WSSR)/WSSR   : 0 delta(WSSR) : 0                 limit for stopping : 1e-05 lambda    : 1150.73  initial set of free parameter values  ...  After 17 iterations the fit converged. final sum of squares of residuals : 3341.93 rel. change during last iteration : -5.29173e-06  degrees of freedom (ndf) : 47 rms of residuals      (stdfit) = sqrt(WSSR/ndf)      : 8.43237 variance of residuals (reduced chisquare) = WSSR/ndf : 71.1049  Final set of parameters            Asymptotic Standard Error =======================            ==========================  a               = 0.26191          +/- 0.005759     (2.199%) b               = 0.00251445       +/- 0.0008358    (33.24%) c               = 0.00541346       +/- 0.0009206    (17.01%) d               = 0.182469         +/- 0.007329     (4.016%)   correlation matrix of the fit parameters:  a      b      c      d a               1.000 b               0.042  1.000 c              -0.229  0.783  1.000 d               0.210 -0.538 -0.768  1.000 gnuplot> replot 
plotfunc3

このように,Yの値が小さい部分以外はフィッティングがあまり良くありません.
これは定義したフィッティング関数が,このデータに対してあまり良くなかったことに
起因します.ピークの形状はLorentzianで良さそうですので, d/sqrt(x) の部分
を,新しいパラメータ e を使って d*x**e のように修正し,その最適値を
gnuplotで求めてみましょう.eの初期値は-0.5としておきます.

 gnuplot> e=-0.5 gnuplot> f(x)=c/((x-a)*(x-a)+b)+d*x**e gnuplot> fit f(x) "exp.dat" using 1:2:3 via a,b,c,d,e  ... Final set of parameters            Asymptotic Standard Error =======================            ==========================  a               = 0.25029          +/- 0.002106     (0.8412%) b               = 0.00197707       +/- 0.0002747    (13.89%) c               = 0.00550098       +/- 0.0003662    (6.657%) d               = 0.21537          +/- 0.003743     (1.738%) e               = -0.358371        +/- 0.0115       (3.208%)   correlation matrix of the fit parameters:  a      b      c      d      e a               1.000 b               0.021  1.000 c              -0.078  0.788  1.000 d              -0.110 -0.384 -0.500  1.000 e              -0.304  0.198  0.335  0.381  1.000 gnuplot> replot 

plotfunc4

Best fitまではあと一歩といった感じですね.a=0.24くらいに固定して
おいて, via b,c,d,eのように a 以外をサーチすると,見た目はもう少し
良くなりますが,χ^2は大きくなってしまいます.


up

Postscriptで出力する

実験データ編のように,出来上がった図を
Postscriptにします.実験データには実線で7番のシンボル○を割り当て,関数には
太めの点線を割り当てます.

 gnuplot> set linestyle 1 lt 1 pt 7 gnuplot> set linestyle 2 lt 2 lw 3 gnuplot> set size 0.6,0.6 gnuplot> set term postscript eps enhanced color Terminal type set to 'postscript' Options are 'eps enhanced color dashed defaultplex "Helvetica" 14' gnuplot> set output "exp.ps" gnuplot> plot"exp.dat" using 1:2:3 title "experiment" with yerrors ls 1, >        f(x) title "Lorentzan" with line ls 2 

[3.8/4.0] Ver.3.8以降のgnuplotの場合は,線種を以下のように定義します.

 gnuplot> set style line 1 lt 1 pt 7 gnuplot> set style line 2 lt 2 lw 3 
plotfunc5


up