Pythonでピークフィッテイング
こんにちは!CheMLです.
今回はPythonを使ったピークフィッテイングについての紹介です.
ピークフィッテイングができるソフトやプログラムはたくさんありますが,もちろんpythonでも可能です.Excelのソルバーでも簡単に実装できますが,速度や使いやすさの観点ではpythonの方が断然おすすめです!それでは早速フィッティングを行っていきましょう.
Gaussianフィッティング
今回はXPS(X線光電子分光)のスペクトル(酸素1s軌道)を用意しました.scipy(※1)のoptimizeパッケージに入っているcurve_fit()という関数を使ってGaussianフィッティングを行います.使い方は簡単です.curve_fitの中身に必要な情報を引数として与えるだけでフィッティングができちゃいます.
curve_fit(関数, 入力データx, 入力データy, [初期値])
ちなみにcurve_fit関数の最適化アルゴリズムにはlevenberg-marquardt法というものが使われています.これはGauss-Newton法と最急降下法を組み合わせた手法で,現在幅広く用いられている最適化の手法です.levenberg-marquardt法については今後詳しく解説する予定です.
curve_fitは戻り値としてフィッティングされたパラメータと係数の共分散行列を返します.そのため,空の変数はp_optとp_covの2つ用意しています.あとは結果のプロットと決定係数の計算,そして結果の出力(csv)です.決定係数の計算は手動でやっていますが面倒な方はpandasのcov()メソッドやscikit-learnのmetrics.r2_score()を使ってもいいと思います.
以下がソースコードです.
gist433a98bfdc7f57973afc3295a1f69c39
以下が実行結果です.savefig()で得られたpng画像を載せました.きれいにフィッテイングできていますね.
(ちなみに高解像度な画像を得るにはオプションにdpi=300を追加するといいです.300くらいで論文・レポートには十分だと思います.Transparent=Trueとすることで背景が透明になります.bbox_inches="tight"は画像サイズをグラフサイズに合わせるためのオプションです.)
各種パラメーターは以下にように表示されました。決定係数は0.9962と良好にフィッテイングできていることが分かります。
gista8299f8fa1a8d33d2ff9b6a09fdd4a5f
擬Voigtフィッティング
以上ではpythonでGaussianフィッティングする方法を紹介しました.しかし実際にはXPSのスペクトルはハイゼンベルグの不確定性幅(寿命に関する線幅,τは寿命,ℏはディラック定数)
を持つためGaussianとLorentzianを畳み込んだVoigt関数を使用すべきですよね.(*Voigt関数は収束性は悪いため畳み込みの代わりにGaussianとLorentzianを足し合わせた擬Voigt関数を用いることが多いです.)
上記のソースコードを以下のように改造すれば擬Voigt関数でのフィッティングが行えます.
gist2550909c75de0eee8eb988d8674e579f
以下が実行結果です.Gaussianフィッティングとほとんど差異がありません.
gist9dde479aa46f0c2a03365036129ba2d6
Gauss関数の割合が1.00にフィッティングされていて,全くLorentzianの寄与がないことが分かります.なおGaussianの割合は使用した装置によってかなり異なるので注意が必要です.
まとめ
以上,今回はscipyのcurve_fitを用いたピークフィッテイングでした!
化学実験のデータ解析では頻繁に使える手法ですので,ぜひ有効活用してみてください.
ご意見・ご感想はコメント欄まで.フォローしていただけるとやる気が出ます!
実行環境
Spyder(anaconda3)
python 3.7.9
windows10
*2023年3月現在では最新版のpython3.8.5に対応していないライブラリも多いため,計算化学ユースにおいてはpython3.7の使用をお勧めします。
※1 Scipy: 科学技術計算用pythonライブラリ(今回のようなフィッティングだけでなくスムージングや関数の微分積分など幅広い演算が可能です。本noteでも多用します。)
★筆者のつぶやき★
「XPSスペクトル」はスペクトルを2回言ってるので「XPスペクトル」が正しいけど違和感ありますね...笑