がれすたさんのDIY日記

電子回路、Python、組み込みシステム開発、自作エフェクターを語るblog

Python:lmfitを使ったカーブフィッティングその2

はじめに

今回は前回の続きでlmfitを使ったフィッティングを試して行きたいと思う。
ほとんどここのパクリだがバグの修正含め解説も入れているので勘弁してほしい

Peak fitting XRD data with Python - Chris Ostrouchov

さっそくやって行こうと思うが
実用的なもので使わないと意味がないので今回はXRD(X-ray Diffraction)のピークフィッティングを例にしてフィッティングしてみたいと思う。
X-ray Diffractionは結晶性物質の解析に用いられる分析で結晶性物質の同定や結晶化度などを求めるときに使う

以下に例としてゼオライトのXRDピークを示す。
f:id:gsmcustomeffects:20181231014222p:plain

XRDの解析ではこれらのピークをフィットしたりして使うことが多いと思う
全部をフィッティングして積分したりすると埒が明かないので特徴的な部分を取り出して局所フィッティングを行うことになる。

これを行うことで非晶質ピークと結晶質ピークの面積がわかるのでそれらの割合から結晶化度を求めることが可能である。

このデータの場合は15度~21度らへんのピークを取り出してやると綺麗にできそうなのでこの範囲でのピークをフィッティングするスキームで話を進めていく
f:id:gsmcustomeffects:20181231014000p:plain


次にフィットモデルについてだがPeak-like modelsってのがLMFITには内蔵されている
XRDはピーク波形の集まりなのでそれを利用すると楽にフィットできそうである。
その中でもよく使うものを紹介しておく

GaussianModel


f(x; A, \mu, \sigma) = \frac{A}{\sigma\sqrt{2\pi}} e^{[ {-{(x-\mu)^2}/{{2\sigma}^2}}] }

LorentzianModel


f(x; A, \mu, \sigma) = \frac{A}{\pi} \big [\frac{\sigma}{(x - \mu)^2 + \sigma^2}\big ]

VoigtModel


f(x; A, \mu, \sigma, \gamma) = \frac{A \textrm{Re}/[ w(z)] }{\sigma\sqrt{2 \pi}}
where

\begin{eqnarray*}
    z &=& \frac{x-\mu +i\gamma}{\sigma\sqrt{2}} \\
    w(z) &=& e^{-z^2}{\operatorname{erfc}}(-iz)
\end{eqnarray*}


lmfitのモデル作成

まずはモデルを作る
def generate_model(spec)という関数を定義する
長いのでColabのリンクを貼っておく(リンク先のgenerate_modelって関数)
https://colab.research.google.com/drive/1At5-CseX7jhqPxjSyTk9gbIxC_dFRRU3#scrollTo=EdlAR8AoLbpa&line=26&uniqifier=1

途中でset_param_hintメソッドという新しい項目が出てきます。
これはモデルの最適化を行う際にある程度のヒントを与えるためのメソッドです。
各種説明は以下の通りです。
f:id:gsmcustomeffects:20181231034341p:plain

これ以外はPythonのビルトインメソッドの多様なのでゆっくり読み進めていけば困ることはないと思います。

次にこの関数に与える辞書リストをつくります。
いきなり作るのもあれなので
まずはピーク検出を行うための関数を作る。これにはScipyのfind_peaks_cwtメソッドを利用する。
こちらもColabのリンクで失礼する(update_spec_from_peaks)
https://colab.research.google.com/drive/1At5-CseX7jhqPxjSyTk9gbIxC_dFRRU3#scrollTo=YESpYSpYRyj6&line=2&uniqifier=1

簡単にモデルを与えてあげる

spec = {
    'x': theta[int((15-5.0)/0.02):int((21-5.0)/0.02)],
    'y': peak[int((15-5.0)/0.02):int((21-5.0)/0.02)],
    'model': [
        {
            'type': 'GaussianModel',
            'params': {'center': 26.5, 'height': 10, 'sigma': 1},
        },
        {
            'type': 'GaussianModel',
            'params': {'center': 26.14, 'height': 350, 'sigma': 0.1},
        },

        {
            'type': 'GaussianModel', 
            'params': {'center': 26.8, 'height': 10, 'sigma': 0.08},
        },         
        {
            'type': 'GaussianModel', 
            'params': {'center': 27.12, 'height': 10, 'sigma': 0.1},
        },
        {
            'type': 'GaussianModel',
            'params': {'center': 27.5, 'height': 10, 'sigma': 0.1},
        },
    ]
}

次にピーク検出したものをプロット

peaks_found = update_spec_from_peaks(spec, [0,1,2,3,4], peak_widths=(13,))
fig, ax = plt.subplots()
ax.scatter(spec['x'], spec['y'], s=4)
peaks = np.zeros((len(peaks_found),2))
j = 0
for i in peaks_found:
    
    ax.axvline(x=spec['x'][i], c='black', linestyle='dotted')
    ax.text(spec['x'][i],800, spec['x'][i])
    peaks[j,0] = spec['x'][i]
    peaks[j,1] = spec['y'][i]
    j = j+1

これを使うとある程度のピーク位置が出てくる
f:id:gsmcustomeffects:20190101022423p:plain

そのままフィッティングしてみる

model, params = generate_model(spec)
output = model.fit(spec['y'], params, x=spec['x'])
fig, gridspec = output.plot(data_kws={'markersize':  1})

f:id:gsmcustomeffects:20190101022833p:plain
仮フィッティング

だいたいあっているがフィットしていない部分もある。3,4番目のモデルが当たっていない感があるのでそこのパラメータの与え方を変えてみる

spec['model'][2]['params']['sigma'] = 0.05
spec['model'][3]['params']['sigma'] = 0.05

f:id:gsmcustomeffects:20190101023609p:plain
ベストフィット

綺麗にできたみたいなのでピークを分けてみてみる

fig, ax = plt.subplots()
ax.scatter(spec['x'], spec['y'], s=4)
components = output.eval_components(x=spec['x'])
print(len(spec['model']))
for i, model in enumerate(spec['model']):
    ax.plot(spec['x'], components[f'm{i}_'])

f:id:gsmcustomeffects:20190101023718p:plain
モデルごとのプロット

最後にベストパラメータを出して終わり

  def print_best_values(spec, output):
    model_params = {
        'GaussianModel':   ['amplitude', 'sigma'],
        'LorentzianModel': ['amplitude', 'sigma'],
        'VoigtModel':      ['amplitude', 'sigma', 'gamma']
    }
    best_values = output.best_values
    print('center    model       amplitude      sigma      ')
    for i, model in enumerate(spec['model']):
        prefix = f'm{i}_'
        values = ', '.join(f'{best_values[prefix+param]:8.3f}' for param in model_params[model["type"]])
        print(f'[{best_values[prefix+"center"]:3.3f}] {model["type"]:16}: {values}')
print_best_values(spec, output)

#--------------------------------------------------------------------------------------------
center    model       amplitude      sigma      
[17.460] GaussianModel   :  428.469,    0.129
[20.031] GaussianModel   :  253.663,    0.153
[19.180] GaussianModel   :  210.520,    0.122
[18.621] GaussianModel   :  254.223,    0.148
[15.726] GaussianModel   :  284.501,    0.123

まとめ

lmfitを使ったパラメータアシストフィットを行うことができた。
またモデルに推定ヒントを単一値を与えていたが範囲で与えることも可能なためその例も書いておいた。
今回のコード全容は以下にアップロードしているので皆さんのほうでカスタマイズしてみてほしい
https://colab.research.google.com/drive/1At5-CseX7jhqPxjSyTk9gbIxC_dFRRU3

今度はレーベンバーグマルカート法以外のメソッドで得たモデルパラメータを使った複合モデリングについて書いていきたい