概要
ある、実部と虚部を有する実験データに対して複素関数用いて回帰分析することでフィッティングパラメータを求める。
pythonではscipyというライブラリが（おそらく）よく知られているが、今回は比較的収束しにくい複素関数のフィッティングに対してパラメータがとりうる値の範囲を制限するためにlmfitというライブラリを利用する。
実験ではインピーダンスの測定や複素誘電率などの周波数応答に位相差が生じる場合などで複素関数が現れる。
環境設定
PC: Windows10 (64bit)
python3.6, jupyter notebook,
lmfit,
理論
今回フィッティングに用いる関数として以下の式を使用する。
Z(f) = R0 + \frac{R}{1+ 2\pi i f C R}
以上の式でははR0, C, Rがフィッティングパラメータで各fに対してzは実部と虚部を持つ。iは虚数。
測定などから得られた実験データをこの式を用いてフィッティングする。
lmfitによるフィッティング
まず、lmfitをインストールする（このときlmfitを動かすためにscipyが必要なのでインストールしていない場合はscipyもインストールする）
pip install lmfit
次にjupyter notebookを開いたら、
- 必要なライブラリをインポートする
#-*-coding:utf-8-*-
import numpy as np
import matplotlib.pylab as plt
from scipy.optimize import leastsq
import math
- 実験データの配列numpy配列を生成
f=np.array([999000, 825200, 681600, 562500, 463900, 382800, 316400, 260700, 215800, 177700,146500, 121100, 99610, 82520, 68120])
rez=np.array([30.77, 35.82, 40.44, 44.99, 49.02, 52.40, 55.05, 57.00,58.41, 59.46, 60.12, 60.90, 61.47, 61.88, 62.19])
imz=np.array([-27.42,-27.00,-26.60,-25.07,-23.12,-20.69,-18.26,-15.90,-13.76,-11.89,-10.35,-7.894,-5.138,-3.568,-3.029])
z=rez+1j*imz
※なお実験データをcsvファイルなどから読み込む場合は配列をnumpy配列に変換する必要がある。
f（グラフの横軸になる）、zの実部rez・虚部imzの配列を実験データなどから取得して、z = (実部）＋i(虚部）の形でnumpy配列を生成。
- フィッティング関数を定義する-
def func(f,R0,C,R):
return R0 + R / (1. + 1j*f*2*math.pi*C*R)
- 初期推定値を入力してフィッティングを実行
from lmfit import Model
model=Model(func)
model.make_params(verbose=True)
print('parameter:', model.param_names, 'independent_variables:',model.independent_vars)
model.set_param_hint('R0',min=0,max=10)
model.set_param_hint('C')
model.set_param_hint('R',min=45,max=80)
result=model.fit(z,f=f, R0=5,C=3.5e-9,R=60)
print('-----finished-----')
model.make_params()でパラメータを自動で生成してくれる（この操作でフィッティングパラメータとしてR0,C,Rができる）。model.set_param_hint()で各パラメータの取りうる値を制限できる。min= ~ で最小値、max= ~ で最大値を設定。他にも値を固定することも可能。model.fit()でフィッティングを実行する。このときに初期推定値を入力する。
- 結果を可視化する
print(result.fit_report())
result.plot(datafmt='o',xlabel='frequency/Hz',ylabel='Re Z',parse_complex='real')
result.plot(xlabel='frequency/Hz',ylabel='Im Z',parse_complex='imag')
フィッティング結果とグラフを表示する。グラフはmatplotlibを使用しているらしい。parse_complex=''で実部realか虚部imagを選択する。保存方法が不明なのでwin + Prt Scでスクリーンショットしてグラフを保存している。
まとめ
lmfitを用いて複素関数のカーブフィッティングを行った。カーブフィッティングについては様々な記事が見られたが複素関数について具体例が示されることが少ないようであるので記事を作成した。複素関数のフィッティングでは通常の実関数と比べて注意すべきことが多いので文献・書籍等に目を通しておくことを推奨する。