gnuplotでガウシアンフィッティングしてみる
gnuplotの使い方が分からないので、兄に相談しながら、試行錯誤したメモを残すべく、このブログを開設した。
(Markdownとは何か、その使い方もよくわからないが、兄に手伝ってもらいつつここに記す)
内容はこんな感じ。
data.txt
0.3 11.56766667 0.45 14.35946667 0.6 2.139466667 0.75 18.39 0.9 -2.0082 1.05 14.37173333 (略) 38.1 1.2624 38.25 -0.508733333 38.4 -0.4528 38.55
このデータを gnuplot を使ってグラフ化する。 この記事を参考にさせていただいた。
gnuplotでガウシアン・フィッティング - 天文学的研究メモ
gnuplotのコマンドラインで下記のように打ち込む。
(フィッティング範囲
とは何のか、どう決めていいのか分からないので、いったん、適当に [-40:40] としてみた)
gnuplot> f(x) = a*exp(-(x-b)**2/(2*c**2)) gnuplot> fit [-40 : 40] f(x) "data.txt" using 1:2 via a,b,c iter chisq delta/lim lambda a b c 0 9.5629421333e+06 0.00e+00 1.56e-01 1.000000e+00 1.000000e+00 1.000000e+00 1 9.4806011957e+06 -8.69e+02 1.56e-02 6.368909e+00 -4.120927e+00 1.530827e+01 * 1.4489883553e+07 3.46e+04 1.56e-01 -8.115759e+01 -1.055083e+01 6.573092e+02 (略) 92 6.5831715808e+05 -9.59e+00 1.56e-07 6.182807e+02 1.912471e+01 1.971323e+00 93 6.5830499570e+05 -1.85e+00 1.56e-08 6.177328e+02 1.912465e+01 1.974877e+00 94 6.5830483668e+05 -2.42e-02 1.56e-09 6.179329e+02 1.912467e+01 1.973602e+00 iter chisq delta/lim lambda a b c After 94 iterations the fit converged. final sum of squares of residuals : 658305 rel. change during last iteration : -2.41571e-007 degrees of freedom (FIT_NDF) : 252 rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 51.1109 variance of residuals (reduced chisquare) = WSSR/ndf : 2612.32 Final set of parameters Asymptotic Standard Error ======================= ========================== a = 617.933 +/- 12.96 (2.097%) b = 19.1247 +/- 0.04781 (0.25%) c = 1.9736 +/- 0.04781 (2.422%) correlation matrix of the fit parameters: a b c a 1.000 b 0.006 1.000 c -0.577 -0.005 1.000 gnuplot> stats "data.txt" * FILE: Records: 255 Out of range: 0 Invalid: 0 Blank: 0 Data Blocks: 1 * COLUMNS: Mean: 19.3500 78.6260 Std Dev: 11.0417 176.9760 Sample StdDev: 11.0635 177.3240 Skewness: 0.0000 2.0501 Kurtosis: 1.8000 5.3076 Avg Dev: 9.5624 123.4886 Sum: 4934.2500 20049.6362 Sum Sq.: 126567.3375 9.56315e+006 Mean Err.: 0.6915 11.0827 Std Dev Err.: 0.4889 7.8366 Skewness Err.: 0.1534 0.1534 Kurtosis Err.: 0.3068 0.3068 Minimum: 0.3000 [ 0] -5.0385 [ 42] Maximum: 38.4000 [254] 558.3840 [114] Quartile: 9.7500 -0.0304 Median: 19.3500 1.5495 Quartile: 28.9500 14.5495 Linear Model: y = -0.4124 x + 86.61 Slope: -0.4124 +- 1.007 Intercept: 86.61 +- 22.44 Correlation: r = -0.02573 Sum xy: 3.751e+005 gnuplot> plot f(x), "data.txt"
出力結果
参考記事から、ガウス関数は
f(x) = a*exp(-(x-b)**2/(2*c**2))
という式であることがわかった。
この結果から半値全幅(FWHM)を求めてみる。求め方については、
を参考にさせていただいた。
gnuplotで使ったガウス関数をLaTeX記法で書いてみる。
はてなブログに数式を書く方法については、
を参考にさせていただいた。
から拝借した。数式上で右クリックすると、LaTeX記法のコードが表示されるのでそれをコピーし、 ^2
などは \^2
とエスケープした。
数式のコードは下記のようになっています。
[tex: \displaystyle f(x)=\frac{1}{\sqrt{2\pi\sigma\^{2}}} \exp\left(-\frac{(x-\mu)\^2}{2\sigma\^2} \right) ]
この式と、最初の式を対照すると、σがcであるので、
であることから、 FWHMの値は
gnuplot> print 2 * c * sqrt(2* log(2)) 4.64747788905307
と求まった。
課題は、 フィッティング範囲
の決め方だ。これがよくわからないので、後日調べることにする。
Excelのマクロ (VBA) のプログラミングの仕方なども覚えたい。。。
その他参考記事