- not so Frequently Asked Questions -
update 2004/8/31
|
|
- not so Frequently Asked Questions - update 2004/8/31
|
gnuplot 入門 --- 関数表示編関数の定義式で書けるような関数であれば，gnuplotはその式を計算して表示すること ができます．「式で書ける」というのは，例えばf(x)=a*x+b のような形で 書けることを表し，微分や積分などを含むものはgnuplotだけでは計算でき ません．gnuplotは三角関数やBessel関数，Gamma関数等，多くの関数を内 部に持っており，それら自身，あるいはそれらを組み合わせたグラフを表 示することができます． 関数のグラフそのものを利用することは(数学の授業目的以外には)少ないと 思います(何かあるかもしれません．あったら教えてください)．gnuplotの 関数機能の威力が最も発揮されるのは，データへの関数当てはめ(データ フィッティング)，つまり最小自乗法です．実験データを，直線等の単純な 解析曲線で近似することがしばしばありますが，gnuplotを使えば非線形関 数の当てはめも簡単に行うことができます．この章では関数の定義の方法 と，その関数に含まれるパラメータに対する最小自乗法を説明します． 次式で表される，Lorentz型関数と 1/sqrt(x) の和の関数を表示してみます． 式に含まれるa,b,c,dの4つの変数は，実験データから求められるものとし ます．  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)   関数値を求める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 最小自乗法でパラメータを求めるパラメータ 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 
上で，パラメータ 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
 このように，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
Best fitまではあと一歩といった感じですね．a=0.24くらいに固定して おいて， via b,c,d,eのように a 以外をサーチすると，見た目はもう少し 良くなりますが，χ^2は大きくなってしまいます． 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  |