gnuplotでガウシアンフィッティングしてみる

gnuplotの使い方が分からないので、兄に相談しながら、試行錯誤したメモを残すべく、このブログを開設した。

(Markdownとは何か、その使い方もよくわからないが、兄に手伝ってもらいつつここに記す)

MTF測定用のタングステン版のエッジデータを貰った。

f:id:wanntarou272:20170521035910p:plain

内容はこんな感じ。

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:id:wanntarou272:20170521035648p:plain

参考記事から、ガウス関数

 f(x) = a*exp(-(x-b)**2/(2*c**2)) 

という式であることがわかった。

この結果から半値全幅(FWHM)を求めてみる。求め方については、

半値全幅(FWHM)と分散の関係 [物理のかぎしっぽ]

を参考にさせていただいた。

gnuplotで使ったガウス関数LaTeX記法で書いてみる。

はてなブログに数式を書く方法については、

cocodrips.hateblo.jp

を参考にさせていただいた。


\displaystyle f(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2} \right)

LaTeX記法でのガウス関数は、

sucrose.hatenablog.com

から拝借した。数式上で右クリックすると、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 = 2\sigma\sqrt{2\log{2}}

であることから、 FWHMの値は

gnuplot> print 2 * c * sqrt(2* log(2))                       
4.64747788905307

と求まった。

課題は、 フィッティング範囲 の決め方だ。これがよくわからないので、後日調べることにする。

Excelのマクロ (VBA) のプログラミングの仕方なども覚えたい。。。

その他参考記事