がれすたさんのDIY日記

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

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

実験データのフィッティングについて頻繁に使う機会があったので自分メモとしてまとめておきます。
フィッティングを行うにあたり、Numpy , Scipyには便利なライブラリがあります。

  • Numpy :polyfit
  • Scipy:optimize.leastsq
  • Scipy:optimize.curve_fit

以上3つが結構有名です。簡単なフィッティングならScipy:optimize.curve_fitを使うのが一番楽だと思います。
一例を示すと

Scipy:optimize.curve_fitを使った例

import numpy as np
import scipy.optimize
import matplotlib.pylab as plt

x = np.linspace(0, 2*np.pi, 100)
y1 = np.sin(x)
y2 = 0.5*np.sin(2*x)
y_noise =(y1 + y2)+ 0.08 * np.random.randn(100)

plt.scatter(x,y_noise)

parameter_initial = np.array([0.5, 1.0]) #a, b, c
# function to fit
def func(x, a, b):
    return np.sin(x) +  a*np.sin(b*x)

paramater_optimal, covariance = scipy.optimize.curve_fit(func, x, y_noise, p0=parameter_initial)
print(paramater_optimal)
fitting = func(x,paramater_optimal[0],paramater_optimal[1])

plt.plot(x,fitting,color="red")

f:id:gsmcustomeffects:20181230003006p:plain
ノイズ入り正弦波合成関数のフィッテング

このように数行でフィットまでできるので
データがそろっていてモデルがある程度絞れているならこの方法が一番早いです。

ですがちょっと複雑なことしようとするとlmfitに軍配が上がるかもしれません。

lmfitって何ぞ?

Non-Linear Least-Squares Minimization and Curve-Fitting for Pythonってサブタイトルがついてる通り非線形最小二乗法を用いたモデルフィットのためのライブラリです。

scipy.optimizeの多くの最適化方法を基にして拡張したものでLevenberg-Marquardt methodの移植からプロジェクトが始まり様々な追加がなされた。
列挙すると以下のとおりである

  • Parameter objectsを導入しモデルにアクセスしやすくした
  • フィッティングアルゴリズムの変更を容易に
  • scipy.optimize.leastsqの信頼区間推定の改善
  • Modelクラスとのカーブフィットが改善(scipy.optimize.curve_fitの機能が拡張)
  • 一般的な線形の多くの組み込みモデルが含まれており、すぐに使用可能。

上記に示すとおりにscipy.optimize.curve_fitの上位互換みたいなものである。

lmfit使用感的な奴

  1. 変数の順序について忘れて意味のある命名でパラメータを参照できる(Pythonの辞書機能でアクセス
  2. パラメータの境界の設定が容易
  3. 目的関数を書き直すことなくパラメータを修正可能
  4. パラメータに代数制約をかけれる(拘束

以下のコードを参考に見てみることとする。

from numpy import exp, sin

from lmfit import minimize, Parameters


def residual(params, x, data, eps_data):
    amp = params['amp']
    phaseshift = params['phase']
    freq = params['frequency']
    decay = params['decay']

    model = amp * sin(x*freq + phaseshift) * exp(-x*x*decay)

    return (data-model) / eps_data


params = Parameters()
params.add('amp', value=10)
params.add('decay', value=0.007)
params.add('phase', value=0.2)
params.add('frequency', value=3.0)

out = minimize(residual, params, args=(x, data, eps_data))

ここで注目すべきはParameter オブジェクトである。
このように他の要素を追加することもできる

params = Parameters()
params.add('amp', value=10, vary=False)
params.add('decay', value=0.007, min=0.0)
params.add('phase', value=0.2)
params.add('frequency', value=3.0, max=10)

vary = Falseの場合、フィット時に値が変化するのを防ぎ、min = 0.0の場合、そのパラメーターの値の下限が設定される。 また、作成後に対応する属性を設定して後で実行することもできる。

以下がPythonの強み。辞書アクセスができる

params['amp'].vary = False
params['decay'].min = 0.10

んで最終的にまとめると

  • paramsオブジェクトをコピーして変更することで、モデルとフィット処理に多くのユーザーレベルの変更を加えることが可能
  • データのモデル化に関する情報の大部分は目的関数に取り込まれるが、ユーザーの外部制御が可能となる。
  • 目的関数の作成者ではなく、フィットを実行するユーザーによる制御が可能。


要はモデル作成者じゃなくても脳死でパラメータ制約部分いじいじするだけでそれっぽいフィットが可能ってわけ
普段から何も考えてない僕にはうれしいツールですね


lmfitを使ってフィットしてみる

さっそくlmfitを利用してフィッティングをやって行きましょう。

まずlmfitを導入します。

pip install lmfit

でいけます。
ほかのパッケージ管理を使ってる場合はご自由にどうぞ

その他必要なものは適宜インストールしてください(Google colabの場合はlmfitのみでOK

import os
import math
import random

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import optimize, signal

from lmfit import models

#gaussianを定義
def gaussian(x, A, μ, σ):
    return A / (σ * math.sqrt(2 * math.pi)) * np.exp(-(x-μ)**2 / (2*σ**2))
  
'''fitするためのノイズ入りデータの作成'''
g_0 = [250.0, 4.0, 5.0]
g_1 = [20.0, -5.0, 1.0]
n = 150
x = np.linspace(-10, 10, n)
y = (gaussian(x, *g_0) + gaussian(x, *g_1)) + 0.05*np.random.randn(n)

fig, ax = plt.subplots()
ax.scatter(x, y, s=20)
'''-----------------------------------------------------------'''

'''lmfitの内蔵モデルの使用'''
model_1 = models.GaussianModel(prefix='m1_')#内包モデルのガウシアンを使用
model_2 = models.GaussianModel(prefix='m2_')#内包モデルのガウシアンを使用
model = model_1 + model_2#コンポジットモデルの作成


'''lmfitの内蔵モデルの初期パラメータセット'''
params_1 = model_1.make_params(center=3.0, sigma=1.0)#ガウシアンのパラメータセット
print(params_1)
params_2 = model_2.make_params(center=-7.0, sigma=1.0)#ガウシアンのパラメータセット
print(params_2)
params = params_1.update(params_2)#update methodにより辞書結合
'''
update methodはPythonのビルトインメソッド(もとから入ってる標準関数みたいなもの
param_1とparam_2の辞書をまとめる機能がある
一般的に同名の辞書があるとその中身が上書きされてしまうがprefix='m1_'を付けているため上書きを避けることができる。

'''
output = model.fit(y, params, x=x)
fig, gridspec = output.plot(data_kws={'markersize': 5})

f:id:gsmcustomeffects:20181230195637p:plain
fittingするための元のデータ

フィッティングあとがこのようになります。

f:id:gsmcustomeffects:20181230195715p:plain
fitting後のデータ

フィッティング後の残差についても出してくれるのでうれしいですね。

fit後のベストパラメータは以下で確認可能です

print(output.best_values)

'''-------------------------------------'''
{'m2_sigma': 4.995352955273267, 'm2_center': 4.000094565660837, 'm2_amplitude': 249.88570423583715, 'm1_sigma': 1.0041137260302269, 'm1_center': -4.998780785999739, 'm1_amplitude': 20.079858838765603}

ざっと書きましたがこんなもんです。
モデル指定してフィットパラメータをいじいじするだけで簡単なものはできます。
応用例とかパラメータ拘束は次回書きます。

この辺勉強するときに使うキーワード

勉強するときに何からやっていいのって方に向けてこの辺で使われるキーワードあげておきます。