がれすたさんの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使用感的な奴

次に使用感だが本家説明が結構しっかりしているので強みみたいのはそれを拝借

まず他のライブラリの問題点をあげる

  • Scipy:optimize.leastsのようにユーザーが引数の順番を理解する必要がある
  • ユーザーがパラメータの拘束をしたい場合に残差の式をいじって対応する必要がある。(できなくはないけどやりにくい
  • パラメーター範囲の設定の時にモデルごとにそれを変更していくのは面倒であり可視的に不利になりやすい

この問題は配列を扱う元実装であるFortranの特性のままであるがPythonのようなリッチなオブジェクトを使うものにはあまり向いていない
そこでパラメーターオブジェクトを扱うことで以下の恩恵が得られる

  1. 変数の順序について忘れて意味のある命名でパラメータを参照できる
  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}

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

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

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