この投稿ではスペクトルなどのデータのカーブフィッティングを、pythonのlmfitというパッケージで行う方法について簡単に書きたい。 ついでなので復習も兼ねて、カーブフィッティングのバックグラウンドについても書いて見た。 lmfitについての詳細は公式ホームページが詳しいのでそちらもご覧いただきたい。
カーブフィッテイングのバックグラウンド
良い機会なのでカーブフィッティングのバックグラウンドについても少しまとめて見ることにした。 数式を使わずに書くので概略になる。 詳しい内容を知りたい人はウェブなり成書を当たって欲しい。 間違った内容・表現を見つけた方がいらっしゃったら、是非コメントでご教授いただきたい。 またlmfitのスクリプトのサンプルを見に来たという方は、投稿下部までジャンプしていただけたら。
さてカーブフィッティングのアイディアはとてもシンプルだ。 あるデータのプロットを適当な曲線で近似するだけである。 しかしその適当な曲線をどう得るかというと、そんなに簡単でもない。
曲線は関数や関数の和として考えることができる。 しかしどう関数を選択するかは自由なので、関数の取り方は無限に存在する。 基本的にはスペクトルなどデータの由来や、カーブの形からシンプルな関数を想定してフィットしていくことになるが、良くフィットしたからといって選んだ関数が正しいとは限らない。
想定する関数が決まったら、その関数を使って何かしらのデータをフィットするわけだ。 関数の種類はすでに想定済みなので、そのパラメータを決めていくことになる。
そのためによく使われるのが最小二乗法だ。 最小二乗法は各データポイントでの関数と実データの誤差の二乗を計算し、その総和が最小になるように関数のパラメータを決める方法だ。
曲線が単純な一つの多項式で表せる場合など関数がパラメータに対して線型結合である場合（関数が線形関数である必要はない）、通常の最小二乗法で正規方程式を解くことでパラメータが得られる。
一方でパラメータが二乗だったり積で表されるような線型結合でない場合は、非線形最小二乗法を使うことになる。 この計算は反復解法と呼ばれるものであり、私はLevenberg-Marquardt法を適用している。 一番よく使われている手法じゃないだろうか。
こうして得られる非線形最小二乗法の解は一つに定まらないので、もっともらしい結果を得るには初期パラメータの設定も重要になる。 スペクトルなどのピークをフィットする場合、ほとんどの場合で非線形最小二乗法を用いてフィットしていくことになる。
カーブフィッティングのためのpythonプログラミング
pythonで非線形最小二乗法のプログラムを書きたい場合だが、scipyに関数が存在している。 Scipyは最も有名というか一般的なpythonの科学用のライブラリである。 Scipyのoptimizeという関数の中に、leastsqとcurve_fitという関数が用意されている。 これらの関数を利用すれば、pythonとしては割と高速のカーブフィットができるかと思う。
例えば無機結晶の粉末回折パターンの解析がしたいというならばこれで十分良いかと思う。 若干関数を定義したりするのがめんどくさかったりもするが大した手間ではない。 まあ粉末解析は各種優秀なソフトが大量にあるので、自分でスクリプトを書く必要なんてのは全くないのだが。
さてなんでlmfitを使うのかという話だが、低結晶性の有機結晶を扱う上でいくつかのlmfitの機能が必要だったのだ。 一番大きな理由はフィッティングパラメータに制限をかけたかったこと。 高結晶性の素材を使ってる人たちに言わせれば、パラメータに制限をかけなきゃいけないなんてのは何かが間違っている証拠てことだそうだ。 しかし低結晶性の素材からそれなりにリーズナブルな値を得るというのも、現実問題として結構需要のあるプロセスなのである。
lmfitを使ったスクリプトの書き方
さてここからは一例を。 まずはlmfitのインストールをpipから。 ちなみに私の環境はanacondaのpython3（anaconda4.4.0 python3.6.1）。
pip install lmfit
実際のコードではまずlmfitのModelを呼び出しておく。 他の必要な関数も適時インポートするがここでは省略。
from lmfit import Model
私の場合はメイン関数の中でデータの読み込みと、データのフィット関数の呼び出し、フィッティング結果を描画する関数の呼び出しをする。 そしてデータのフィット関数と描画の関数を定義するが、ここら辺は好みで。 データは何かしら２カラムで、x軸に対して強度を取ったようなものを読み込むことが多いかな。
plot = np.loadtxt(argvs[1]) comps, result, int, cent = Fitex(plot) namefig = argvs[1] +".png" plot_fitted(plot,namefig,comps,result)
続いてフィッティングに使う関数の定義。 この例では面積計算のガウシアンとコンスタントを定義。
def gaussian(x, amp, cen, hf): return amp*(2/hf)*np.sqrt((np.log(2)/np.pi))*np.exp((-4*np.log(2))*((cen-x)/hf)**2) def consta(x,con): return con
フィッティング操作を行う関数は以下のような感じ。
def Fitex(arr):
＃x軸とy軸のカラムの指定。 もちろん自分のデータ次第。
y = arr[:,1] x = arr[:,0]
＃個別のフィッティング関数の設定。 ここではガウシアンとコンスタント。 私はフィッティング関数は自分で定義しているけど、lmfitはプリインストールモデルが結構あるはず。 prefixは後でパラメータ設定する用の名前付け。
gauss1 = Model(gaussian,prefix='g1_') Const1 = Model(consta,prefix='c1_')
＃フィットする関数群を決める。 ガウシアンとコンスタントの和。
mod = gauss1 +Const1
＃フィッティングに使うパラメータを定義
param = gauss1.make_params() param.update(Const1.make_params())
＃フィッティングパラメータの初期値や制限などの設定。 ここを結構自由に設定できるのがlmfitの良いところ。 例えば、ここではConst1関数は０から１３００の範囲での最小値を拾って初期値とし、その値はフィッティングの間不変に設定されている（pick_minは自己定義関数）。 gauss1関数はフィッティングに使われるが、min,maxで最大値・最小値が指定されている。 あとはexprを使って式で値を決定することもできる。 例えばparam[‘c2_con’].set(expr=’c1_con*3′)といった具合。
param['c1_con'].set(pick_min(y[0:1300]),vary=False) param['g1_amp'].set(100000, min=0, max=5000000) param['g1_cen'].set(180,min=170, max=190) param['g1_hf'].set(200,min=1, max=3000)
＃フィッティングの実行。 結果はresultに収納。
result = mod.fit(y, param, x=x)
＃あとは適当にフィッテイング結果の出力
print(result.fit_report()) print(result.params)
＃フィッティングの結果の値を他の関数から呼び出したりして使いたい時は、best_values.getで取得しておく。
comps = result.eval_components(x=x)
int = result.best_values.get('g1_amp')
cent = result.best_values.get('g1_cen')
＃最後に使いたい値をリターンする。
return comps,result, int, cent
以上lmfitを使ったカーブフィッティングの例でした。 後はresult.best_fitや出力のcomps[‘g1_’]+comps[‘c1_’]などを適当なx軸に対してプロットしてあげれば絵も描ける。
最後にlmfitの公式ホームページのアドレスと、使用した際の引用先を記しておく。
[1]https://lmfit.github.io/lmfit-py/
[2]http://dx.doi.org/10.5281/zenodo.11813
関連記事
D