Gnuplotでカーブフィッティング

実験で得られた数値スペクトルデータがあるとする。
これには、2つのガウシアン型のピークを持つ成分が含まれている。
ガウシアン型は $ y = A_1 \times \exp{[A_2 (x-A_3)^2]} $ 型。
この実験結果から成分を抽出する方法を解説します。

Gnuplotでのフィッティングコマンドは

# 関数定義
f(x) = a*x + b  # 任意の関数・変数
    a = 1.0     # 初期値
    b = -5.0    # 

# フィッティング
fit [xmin:xmax] f(x) "解析データ" u 1:2 via a, b

Step. 1 とりあえずプロットする

fort.100 というのが実験データ( 1列目 : x, 3列目: y )とする。

plot "fort.100" u 1:3 

見ると、バックグラウンドと $x = 8, 13$あたりにピークがありそう。
振幅はどちらも2-4くらいに見える(およその目安把握が重要)。

Step. 2 バックグラウンドを計算して、結果から差し引く

# y = ax + b 型のバックグラウンド処理

B(x) = B1*x + B2
    B1 = 0.3
    B2 = 1.0

fit [17:19] B(x) "fort.100" u 1:3 via B1, B2

plot \
"fort.100" u 1:3 w l lc "blue" t "Original",\
B(x) w l lc "black" dt (5,3) t "Back ground",\
"fort.100" u 1:(($3)-(B1*($1)+B2)) w l lc "red" lw 2 t "subtract back ground"

バックグラウンドは直線に乗ってそうな $x = 17-19$ でフィッティング。
実行結果は以下の通り。
問題なくバックグラウンド処理はできていそう。

Step. 3 ガウシアンフィッティング

# 2つのガウシアン波形を含む h(x) でフィッティング

h(x) = f1*exp( f2*(x-f3)**2 ) + g1*exp( g2*(x-g3)**2 )
    f1 = 1.0        # 振幅 : 目視レベルで入力
    f2 = -1.0       # 横幅 : ガウシアンなので負の値入力
    f3 = 8.0        # 位置 : 目視レベルで入力
    g1 = 1.0        
    g2 = -1.0
    g3 = 13.0

fit [4:15] h(x) "fort.100" u 1:(($3)-(B1*($1)+B2)) via f1, f2, f3, g1, g2, g3

plot \
"fort.100" u 1:(($3)-(B1*($1)+B2)) w l lc "black" lw 8,\
h(x) w l lw 2 lc "yellow" dt (10,2) t "fitted curve",\

結果は次の通り

B1              = 0.0999337
B2              = 2.00122

f1              = 2.99945
f2              = -0.250139
f3              = 8.00014
g1              = 0.749658
g2              = -0.500211
g3              = 13

実際は $y = 3 e^{-0.25(x-8)^2} + 0.75 e^{-0.5(x-13)^2} + 0.1x + 2$ を計算していたので、精度よくフィッティングができたといえる。

注意点

初期値は適切な値を入れること