見出し画像

【Python】lmfitを少しでも使いやすくしてみたかった。

のうない

こんにちは。のうないです。

サイエンスの世界では、観測されたデータを数式に当てはめてパラメータを算出する、ということを頻繁にやります。
俗にいう、カーブフィッティングです。

そのカーブフィッティングをするための、lmfitというライブラリがあります。lmfitはscipy.optimizeライブラリの強化版とも言えるものです。
線形関数やガウス関数をはじめ、多くのモデル関数があらかじめ用意されている他、自分で関数を定義することもできます。

非常に汎用性が高いのですが、ただ、一介の実験系研究者ののうないにとっては、ちょっと取っ付きにくいなと思うところもあります。
基本的な使い方は、以下の手順になります。

  1. モデルオブジェクトを生成する。

  2. パラメータオブジェクトを生成し、初期パラメータをセットする。

  3. フィッティングを実行する。

  4. 結果を取り出す。

一例を示します。
例えば、とあるデータに対して、ガウス関数３つと一次関数の足し算でフィッティングしたいと思ったとします。

1. モデルオブジェクトの生成

from lmfit import models

# モデルオブジェクトを生成
gauss1 = models.GaussianModel(prefix='gauss1_') #ガウス関数1
gauss2 = models.GaussianModel(prefix='gauss2_') #ガウス関数2
gauss3 = models.GaussianModel(prefix='gauss3_') #ガウス関数3
linear = models.LinearModel(prefix='linear_')   #一次関数

# モデルオブジェクトを足す
total_model = gauss1 + gauss2 + gauss3 + linear

2. パラメータオブジェクトの生成と初期パラメータの設定

# 各モデル関数のパラメータオブジェクトを生成
gauss1.param = gauss1.make_params(center=dict(value=100, min=125, max=150),
                                  sigma=dict(value=15, min=0),
                                  amplitude=dict(value=2000, min=0))
gauss2.param = gauss2.make_params(center=dict(value=400, min=300, max=500),
                                  sigma=dict(value=20, min=0, max=35),
                                  amplitude=dict(value=17000, min=0))
gauss3.param = gauss3.make_params(center=dict(value=600, min=550, max=650),
                                  sigma=dict(value=10, min=0, max=15),
                                  amplitude=dict(value=1000, min=0))
linear.param = linear.make_params(slope=dict(value=3, min=0, max=8),
                                  intercept=dict(value=0, min=-10, max=10))
# パラメータオブジェクトを足す
total_param = gauss1.param + gauss2.param + gauss3.param + linear.param
total_param

3.フィッティングを実行する

fit_result = total_model.fit(data=vector_y, params=total_param, x=vector_x)
# vector_x, vector_yはフィッティングしたいX軸、Y軸データ。

4.結果を取り出す、は省略。
モデル関数が複雑だと、どうしてもコード量が多くなってしまいます。
特に、モデル関数の係数名(centerとか、interceptとか)を覚えておくのが個人的にメンドイです。

lmfitを継承して新しいクラスを作ってしまう

そこで、のうないはlmfitをベースに新しくクラスを作ることを考えました。
ここではガウス関数の例を示します。

from lmfit.models import GaussianModel

# lmfit.model.GaussianModelを継承して、新しくクラス定義してみる
class MyGaussian(GaussianModel):
    def __init__(self, prefix:str='', height:float=None, center:float=None, sigma:float=None):
        super().__init__(prefix=prefix) 

        # インスタンス時にパラメータオブジェクトを生成しておく。
        self.params = self.make_params()
        self.params[f'{prefix}amplitude'].set(value=height/0.4*sigma, min=-inf, max=inf, vary=True)
        self.params[f'{prefix}center'].set(value=center, min=center-sigma, max=center+sigma, vary=True)
        self.params[f'{prefix}sigma'].set(value=sigma, min=0, max=sigma*5, vary=True)
        
        # フィッティングの結果を格納するためのフィールド
        self.fitted_y = None
        self.amplitude = None
        self.center = None
        self.sigma = None
        self.area = None
        self.fwhm = None
    
    def hogehoge(): #お好みのメソッドを実装
        pass

新しくクラスを作るのもメンドイですが、メリットは多いです。

  • モデルオブジェクトの生成時にパラメータオブジェクトも同時に生成する

  • モデル関数の係数名を覚えておく必要がない

  • 係数の初期値はモデルオブジェクト生成時の引数で渡す

  • モデル関数独自の追加処理を実装できる(面積計算、パラメータ変換など)

このクラスを使うと、先ほどのコードがこうなります。

# モデルオブジェクトを生成
gauss1 = MyGaussian(prefix='gauss1_', height=xxxx1, center=yyyy1, sigma=zzzz1) #ガウス関数1
gauss2 = MyGaussian(prefix='gauss2_', height=xxxx2, center=yyyy2, sigma=zzzz2) #ガウス関数2
gauss3 = MyGaussian(prefix='gauss3_', height=xxxx3, center=yyyy3, sigma=zzzz3) #ガウス関数3
linear = MyLinear(prefix='linear_', slope=aaaa, intercept=bbbb) #一次関数

# モデルオブジェクトを足す
total_model = gauss1 + gauss2 + gauss3 + linear
# パラメータオブジェクトを足す
total_param = gauss1.params + gauss2.params + gauss3.params + linear.params

# フィッティングの実行
fit_result = total_model.fit(data=vector_y, param=total_params, x=vector_x)

手順自体は変わりませんが、ゴリッゴリのプログラマではない人が普段使いする分には使いやすくなったのではないでしょうか。

いいなと思ったら応援しよう！