scipyを使ったピークサーチ、最適化のやり方について紹介します。

この記事、X線回折やFTIRなど、実験データにフィッティングして数値を求める場合を想定して書いています。

この記事を書いた背景

(興味ない方は、読み飛ばしてもらってもいいです)

ほとんどの場合、実験データへのフィッティングは測定機に付属している専用ソフトウェアを使用すると思います。
しかし、データによっては自動でのフィッティングが難しかったり、バッチ処理ができずデータを一つずつしかフィッティングできない、などの場合があります。

だったら、pythonでフィッティングプログラムを自作すればいいじゃないか！
ということで、自作した際に一番苦労したピークサーチとパラメータ最適化について、まとめておこうと思いました。

scipyでピークサーチ

まずは、scipyでピークサーチするやり方について紹介します。
二つのやり方を紹介します。
①scipy.optimizeのcurve_fitを使ってフィッティングする
②scipy.signalのargrelmaxを使って極大値をみつける
となります。

①はフィッティングする関数が分かっているときに有効で、
②はフィッティング不要でピークをみつけることができます。

scipy curve_fitを使う場合

フィッティングする関数を定義

フィッティングする関数を定義します。
ピークへのフィッティングにはガウス関数やローレンツ関数などが使われます。

ガウス関数

x = bで最大値 aをとります。

ローレンツ関数

x = bで最大値aをとります。また、dは半値全幅(HWHM)になります。

ローレンツ関数であれば、次のように定義します。

def lorentz_func(x, A, X0, HWHM):
return A * HWHM / ((x - X0)**2 + HWHM**2)

フィッティングパラメータはA, X0, HWHMの３つです。

scipy curve_fitを使ってフィッティング

scipy.optimizeにあるcurve_fitを使ってフィッティングします。

curve_fitに、フィッティングする関数、フィッティングしたいデータ(x, y)、フィッティングパラメータの初期値を変数として、代入します。（最低限この4つ）
また、curve_fitは返り値として最終的なフィッティングパラメータと共分散を返します。

コードは次のように書きます。

param = [A_ini, X0_ini, HWHM_ini] #パラメータ初期値
popt, pcov = curve_fit(lorentz_func, om, inte, p0=param)

poptが最適化されたパラメータになります。

例えば、ピーク値をとるxの値が知りたければpopt[1]をみればよいです。

X線回折波形のピークにローレンツ関数をフィッティングした結果。
青色：フィッティング結果、橙色：実験データ

フィッティングする際の注意事項

フィッティングするうえで、パラメータの初期値の選び方が重要です。初期値があまりにも外れすぎているとうまくフィッティングできないので、実験データを見て、適切な初期値を設定しましょう。

scipy argrelmaxを使う場合

scipy.signalの極大値を求める関数、argrelmaxを使うと、
curve_fitよりも簡単です。

signal.argrelmax(y, order=3)

で極大値のインデックスを返してきます。
yは極大値を探したいデータ、orderは極大値を探す間隔です。
例えば、order=1としてしまうと、ノイズが多いデータの場合、非常に多くの点がピークと判断されてしまいます。
ちょうどいい値に調整しましょう。

scipyで実験データへのフィッティング

続いて、実験データへのフィッティングのやり方を紹介します。
すでに紹介したcurve_fitが実験データへのフィッティングに使えます。
ですが、ここではそれとは違う方法を紹介します。

scipy least_squaresを使ったフィッテイング

scipy.optimizeのleast_squaresを使ってフィッティングしてみます。この関数は最小二乗法でフィッティングする関数です（curve_fitも最小二乗法を使っています）。

least_squaresと同じような関数として、leastsqという関数もあります。二つの違いは何かというと、leastsqはパラメータの範囲を制限できないのに対して、least_squqresではパラメータの範囲を指定できるということです。
実験データのフィッティングではある程度パラメータの取りうる範囲が決まっていますので、least_squaresを使用するのが良いと思います。

least_squaresの使い方

まずは、フィッティングする関数を定義します。
curve_fitと違うのは、フィッティングする関数と実験データの差分の関数を定義します。

例えば、

def fitting_func(x, a, b, c):

f = a*x**2 + b*x + c

return f

def ls_res(a, b, c, x, y):

r = y - fitting_func(x, a, b, c)

return r

のように定義します。

次にパラメータの範囲を定義します。

パラメータの範囲は、それぞれの上限と下限を配列で定義します。
例えば、

の場合、

b = ([-1, 0, -np.inf], [1, np.inf, np.inf])

となります。

そして、

param = [a_ini, b_ini, c_ini]

b = ([-1, 0, -np.inf], [1, np.inf, np.inf])

le_lsq = least_squares(les_res, param, args=(x, y), bounds=b)

でフィッチングが開始されます。
フィッティング結果は、

popt = le_lsq.x

で取り出せます。

まとめ

実験データのピークサーチのやり方、フィッティングの方法を紹介しました。
皆さんが自作でフィッティングソフトを作る時の参考になれば幸いです。

参考ページ

ローレンツ関数でのフィッティング
https://qiita.com/yokotani/items/f8920f65b1037ec7009d

argrelmaxでピークの抽出
https://qiita.com/wrblue_mica34/items/e174a71570abb710dcfb

X線回折波形をpythonで計算
https://sciencompass.com/phys-engineer/semiconductor_physics/xrayutilities_1 ‎
https://sciencompass.com/phys-engineer/semiconductor_physics/xrayutilities_2