こんな人にオススメ
pythonを使ってデータを関数とか既存の他のデータにフィットさせたい（フィッティング）。どうしたらいい？
ということで今回はpythonを使ってデータをフィッティングする方法について解説する。pythonでのフィッティングの方法はいくつかあるが、今回はnumpyのpolyfitとscipyのcurvefitにフォーカスする。
python環境は以下。
- Python 3.9.4
- numpy 1.20.3
- scipy 1.6.3
- plotly 4.14.3
- plotly-orca 3.4.2
下準備
import sys
import numpy as np
import scipy.optimize
import plotly.graph_objects as go
import plotly.io as pio
sys.path.append('../../')
import plotly_layout_template as template
まずは今回使用するモジュールのimport。plotly_layout_templateは自作関数で、2つ上のディレクトリに置いているので../../で移動している。詳細は以下参照。
-
-
【随時更新 備忘録】plotlyのグラフ即席作成コード
こんな人にオススメ plotlyっていじれ& ...
続きを見る
1, 2次関数をfit（np.polyfit）
まずはnumpyのpolyfit。これは1次関数や4次関数などの多項式をfitするための関数で、多項式なら簡単にfitすることができる。一方で複雑な関数には対応できないので、その場合は後述のscipy.optimize.curve_fitを使用する。
グラフ作成の準備
def scatter(mode, x, y, name, symbol):
"""plotlyのScatterの雛形関数
Parameters
----------
mode : str
プロットの種類
x : array
横軸の値
y : array
縦軸の値
name : str
プロットの名前
symbol : str
プロットのマーカー
Returns
-------
plotly.graph_objs._scatter.Scatter
plotlyのScatterデータ
"""
d = go.Scatter(
mode=mode,
x=x, y=y,
name=name,
marker=dict(symbol=symbol),
hovertemplate=f"{name}<br>"
+ 'x: %{x} <br>'
+ f"{name}" + ': %{y} <br>'
+ "<extra></extra>",
)
return d
今回はfitさせたのち、元のプロットとfit後のプロットをグラフ化する。そのために関数を2種類用意する。ますはデータのプロットを表すscatter関数。データ配列x, yを入れるだけでplotlyのプロット配列を作成することができる。
また、グラフ自体を作成するためにgraph関数も定義した。
def graph(x, y, fit, title: str, save_name: str):
"""fitしたい、fit後のプロットをグラフ化
Parameters
----------
x : array
横軸の値
y : array
fitしたいデータ
fit : array
fit後のデータ
title : str
グラフタイトル
save_name : str
ファイル保存名
"""
y_dct = {'y': y, 'fit': fit}
m_dct = {'y': 'markers', 'fit': 'lines+markers'}
s_dct = {'y': 'circle', 'fit': 'triangle-up'}
plot = []
for name, val in y_dct.items():
d = scatter(
mode=m_dct[name],
x=x, y=val, name=name,
symbol=s_dct[name],
)
plot.append(d)
layout = go.Layout(
template=template.plotly_layout(),
title=dict(text=f"{title}",),
xaxis=dict(title='x',),
yaxis=dict(title='y',),
)
fig = go.Figure(data=plot, layout=layout)
fig.show(config=template.plotly_config())
pio.write_html(
fig, f"{save_name}.html",
auto_open=False, config=template.plotly_config()
)
pio.write_image(fig, f"{save_name}.png")
y_dctで元のデータとfitしたデータを1つのdictに入れてそれをプロットする。ここでさっきのscatter関数を使用してコードを短めにしている。
データの作成
x = np.arange(20) print(x) # [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19]
横軸の値は上に示すように0から19までを1刻みの合計20データ。縦軸の値は多項式なので、以下のget_polyfit関数でデータを作成、さらにfitまで行う。引数のpowerは何次関数のデータを作成するのか、orderはそれを何次関数でfitするのかを示す。
また、c_tfはfitの共分散行列を、f_tfはランクなどのデータをfit時に取得するかどうかを表す。なお、c_tfとf_tfは同時に使用することができない。
def get_polyfit(power: int, order: int, c_tf=False, f_tf=False):
"""線形関数を引数で作成
作成したデータをnumpyのpolyfitで線形fit
Parameters
----------
power : int
作成データの次数
order : int
fitの次数
c_tf : bool, optional
分散行列を出力するか否か, by default False
f_tf : bool, optional
パラメータ、残差、ランク、
特異点の値？、特異値のカットオフ比？を出力するか否か, by default False
Returns
-------
[type]
tuple
作成データ
fitパラメータ
fit結果
"""
# データの作成
y = []
for i in range(20):
np.random.seed(seed=i) # 各データでseedを変更する
val = i + np.random.rand() # 乱数によるズレを大きくする
y.append(val ** power)
y = np.array(y)
# numpyで線形fitし係数を係数を計算
if c_tf: # パラメータ、共分散を出力
par, cov = np.polyfit(x, y, order, cov=True)
fit = np.poly1d(par)(x) # 係数からfit配列を作成
return y, par, cov, fit
elif f_tf: # パラメータ、残差、ランク、特異点の値？、特異値のカットオフ比？を出力
par, res, rank, singular, rcond = np.polyfit(x, y, order, full=True)
fit = np.poly1d(par)(x)
return y, par, res, rank, singular, rcond, fit
else: # パラメータだけ出力
par = np.polyfit(x, y, order)
fit = np.poly1d(par)(x)
return y, par, fit
多項式データはxごとに乱数を別に作成して一定の乱数値にならないようにしている。また、今回はf_tfは使用しないが、何かしら必要な時は使用するといいだろう。
1次関数っぽいデータを1次関数でfit
# 1次関数っぽいのを1次関数でfit
y, par, fit = get_polyfit(power=1, order=1)
print(f"y: {y}")
# y: [ 0.5488135 1.417022 2.4359949 3.5507979 4.96702984 5.22199317
# 6.89286015 7.07630829 8.8734294 9.01037415 10.77132064 11.18026969
# 12.15416284 13.77770241 14.51394334 15.8488177 16.22329108 17.294665
# 18.65037424 19.0975336 ]
print(f"par: {par}")
# par: [0.99143382 0.55671386]
print(f"fit: {fit}")
# fit: [ 0.55671386 1.54814769 2.53958151 3.53101533 4.52244916 5.51388298
# 6.50531681 7.49675063 8.48818446 9.47961828 10.47105211 11.46248593
# 12.45391975 13.44535358 14.4367874 15.42822123 16.41965505 17.41108888
# 18.4025227 19.39395653]
まずは1次関数っぽいデータ配列を1次関数でfitする。既に関数を定義しているので、引数に1次関数のデータを作成、1次関数でfitという情報を入れると完成。実際にグラフ化する時は以下のコードでグラフ化する。
# グラフ作成
graph(
x=x, y=y, fit=fit,
title='linear', save_name='linear_polyfit',
)
2次関数っぽいデータを1次関数でfit
# 2次関数っぽいのを1次関数でfit
y, par, fit = get_polyfit(power=2, order=1)
print(f"y: {y}")
# y: [3.01196262e-01 2.00795136e+00 5.93407116e+00 1.26081657e+01
# 2.46713854e+01 2.72692127e+01 4.75115211e+01 5.00741390e+01
# 7.87377494e+01 8.11868424e+01 1.16021348e+02 1.24998430e+02
# 1.47723674e+02 1.89825084e+02 2.10654551e+02 2.51185022e+02
# 2.63195173e+02 2.99105438e+02 3.47836459e+02 3.64715790e+02]
print(f"par: {par}")
# par: [ 19.74311333 -55.28141636]
print(f"fit: {fit}")
# fit: [-55.28141636 -35.53830303 -15.7951897 3.94792362 23.69103695
# 43.43415028 63.17726361 82.92037694 102.66349026 122.40660359
# 142.14971692 161.89283025 181.63594357 201.3790569 221.12217023
# 240.86528356 260.60839689 280.35151021 300.09462354 319.83773687]
お次は2次関数っぽいデータだが、初めは1次関数でfitすることとする。2次関数っぽいからといって必ず2次関数でfitする必要はない。時と場合によっては1次関数でfitするのが正しい時もあるかもしれない。グラフ化以下。
graph(
x=x, y=y, fit=fit,
title='linear (quad)', save_name='linear_quad_polyfit',
)
2次関数っぽいデータを2次関数でfit
# 2次関数っぽいのを2次関数でfit
y, par, fit = get_polyfit(power=2, order=2)
print(f"y: {y}")
# y: [3.01196262e-01 2.00795136e+00 5.93407116e+00 1.26081657e+01
# 2.46713854e+01 2.72692127e+01 4.75115211e+01 5.00741390e+01
# 7.87377494e+01 8.11868424e+01 1.16021348e+02 1.24998430e+02
# 1.47723674e+02 1.89825084e+02 2.10654551e+02 2.51185022e+02
# 2.63195173e+02 2.99105438e+02 3.47836459e+02 3.64715790e+02]
print(f"par: {par}")
# par: [ 0.96144676 1.47562492 -0.47895115]
print(f"fit: {fit}")
# fit: [ -0.47895115 1.95812053 6.31808573 12.60094445 20.80669668
# 30.93534243 42.98688169 56.96131447 72.85864076 90.67886058
# 110.4219739 132.08798075 155.67688111 181.18867498 208.62336238
# 237.98094328 269.26141771 302.46478565 337.59104711 374.64020208]
続いて2次関数っぽいデータを2次関数でfitする。先ほどと異なるのはfitの次数を表すorderを2にするというとこ。これに伴ってパラメータは3つになる。グラフは以下。
graph(
x=x, y=y, fit=fit,
title='quad', save_name='quad_polyfit',
)
適当に作成した3次関数を3次関数でfit
# 適当に3次関数
y = 0.05 * x ** 3 - 1.1 * x ** 2 + 10 * x ** 1 + 1.5 * x ** 0
par = np.polyfit(x, y, 3) # numpyで線形fitし係数を係数を計算
fit = np.poly1d(par)(x) # 係数からfit配列を作成
print(f"y: {y}")
# y: [ 1.5 10.45 17.5 22.95 27.1 30.25 32.7 34.75 36.7 38.85
# 41.5 44.95 49.5 55.45 63.1 72.75 84.7 99.25 116.7 137.35]
print(f"par: {par}")
# par: [ 0.05 -1.1 10. 1.5 ]
print(f"fit: {fit}")
# fit: [ 1.5 10.45 17.5 22.95 27.1 30.25 32.7 34.75 36.7 38.85
# 41.5 44.95 49.5 55.45 63.1 72.75 84.7 99.25 116.7 137.35]
最後は適当に作成した3次関数を3次関数でfit。まあ同じ。パラメータは4つ。
graph(
x=x, y=y, fit=fit,
title='cubic', save_name='cubic_polyfit',
)
np.polyfitでパラメータの誤差を求める
fitでパラメータを決めたがこのパラメータには誤差（エラー）がある。ここではnumpy.polyfitでのfitの誤差を求める。
誤差を求めるにはget_polyfit関数の引数c_tfをTrueにして共分散行列を出力。そのあと、そのnp.diagで対角線をとり、その平方根を取ることで誤差が計算できる。共分散行列とか細かいことは省略する。$S_{xx}$と$S_{yy}$を使用している。
1次関数っぽいのを1次関数でfit
y, par, cov, fit = get_polyfit(power=1, order=1, c_tf=True)
print(f"par: {par}")
# par: [0.99143382 0.55671386]
print(f"cov: {cov}")
# cov: [[ 0.00014713 -0.00139776]
# [-0.00139776 0.01817087]]
まずは1次関数っぽいのを1次関数でfitするパターン。共分散行列はcovで表している。これの対角線の平方根を計算してパラメータの誤差を計算可能。
# 共分散行列の対角線の平方根がフィット係数の誤差に該当
err = {}
for num, p in enumerate(np.sqrt(np.diag(cov)), 1):
err[f"err{num}"] = p
print(err)
# {'err1': 0.012129820997481852, 'err2': 0.1347993725605669}
ちなみに本来は有効数字とか考えないといけないが今回は無視。
2次関数っぽいのを1次関数でfit
# 2次関数っぽいのを1次関数でfit
y, par, cov, fit = get_polyfit(power=2, order=1, c_tf=True)
print(f"par: {par}")
# par: [ 19.74311333 -55.28141636]
print(f"cov: {cov}")
# cov: [[ 1.42853734 -13.57110475]
# [-13.57110475 176.42436176]]
# 共分散行列の対角線の平方根がフィット係数の誤差に該当
err = {}
for num, p in enumerate(np.sqrt(np.diag(cov)), 1):
err[f"err{num}"] = p
print(err)
# {'err1': 1.1952143498851289, 'err2': 13.282483267753786}
2次関数っぽいのを1次関数でfitバージョン。
2次関数っぽいのを2次関数でfit
# 2次関数っぽいのを2次関数でfit
y, par, cov, fit = get_polyfit(power=2, order=2, c_tf=True)
print(f"par: {par}")
# par: [ 0.96144676 1.47562492 -0.47895115]
print(f"cov: {cov}")
# cov: [[ 2.91899204e-03 -5.54608488e-02 1.66382546e-01]
# [-5.54608488e-02 1.13081752e+00 -3.89335159e+00]
# [ 1.66382546e-01 -3.89335159e+00 1.90008868e+01]]
# 共分散行列の対角線の平方根がフィット係数の誤差に該当
err = {}
for num, p in enumerate(np.sqrt(np.diag(cov)), 1):
err[f"err{num}"] = p
print(err)
# {'err1': 0.05402769699957549, 'err2': 1.0633990396311224, 'err3': 4.35900066588794}
2次関数っぽいのを2次関数でfitバージョン。2次関数でfitするのでパラメータが3種類になる。これに伴って誤差も3種類になる。
1, 2次関数をfit（scipy.optimize.curve_fit）
上で見たnumpy.polyfitは多項式関数を使用するときにはお手軽でいいんだけど、多項式関数以外の場合はfitできない。まあ名前がpolyfit（polynomial: 多項式）なのでしゃーない。
じゃあ、多項式関数以外はfitできないかというとそうではない。scipy.optimize.curve_fitを使用することで任意の関数にfitさせることが可能。ここからはこのcurve_fitを見ていく。
データ・グラフ作成の準備
def get_poly(power: int):
# データの作成
arr = []
for i in range(20):
np.random.seed(seed=i) # 各データでseedを変更する
val = i + np.random.rand() # 乱数によるズレを大きくする
arr.append(val ** power)
arr = np.array(arr)
return arr
まずはデータ作成の準備としてget_poly関数を定義。polyfitとは異なり、この関数では単純に多項式関数の作成のみ。また、横軸の値は先ほどと同様、0から19まで1刻み。
# 横軸の値 x = np.arange(20) print(x) # [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19]
さらに再度scatter関数、graph関数を使用するので以下に再掲する。確認したい方は開いてほしい。
1次関数っぽいデータを1次関数でfit
def linear(x, a, b):
y = a * x + b
return y
arr = get_poly(power=1)
par, cov = scipy.optimize.curve_fit(f=linear, xdata=x, ydata=arr)
fit = linear(x, *par)
curve_fitを使用するには予めfitしたい関数の形を定義しておかないといけない。今回は1次関数っぽいデータを1次関数でfitするので、1次関数の形状のlinear関数を定義した。
arrで縦軸の値（fitしたいデータ）を作成、scipy.optimize.curve_fitの引数でfitしたい関数と横・縦軸の値を指定すれば良い。戻り値はパラメータと共分散行列。それぞれparとcovに代入した。
あとはlinear関数に横軸のxとパラメータparを展開して入れればfit関数が出来上がる。それぞれの出力値は以下。微妙にpolyfitとズレているので使用する際には桁数に注意が必要。
print(par) # [0.99143383 0.55671385] print(cov) # [[ 0.00014713 -0.00139776] # [-0.00139776 0.01817087]] print(fit) # [ 0.55671385 1.54814768 2.5395815 3.53101533 4.52244915 5.51388298 # 6.5053168 7.49675063 8.48818445 9.47961828 10.4710521 11.46248593 # 12.45391975 13.44535358 14.43678741 15.42822123 16.41965506 17.41108888 # 18.40252271 19.39395653]
グラフ化しても内容は変わらないが、一応グラフ化。
graph(x=x, y=arr, fit=fit, title='linear', save_name='linear_curve_fit')
2次関数っぽいデータを1, 2次関数でfit
arr = get_poly(power=2)
par, cov = scipy.optimize.curve_fit(f=linear, xdata=x, ydata=arr)
fit = linear(x, *par)
print(par)
# [ 19.74311333 -55.28141636]
print(cov)
# [[ 1.4285374 -13.57110512]
# [-13.57110512 176.42436267]]
print(fit)
# [-55.28141636 -35.53830303 -15.7951897 3.94792362 23.69103695
# 43.43415028 63.17726361 82.92037694 102.66349026 122.40660359
# 142.14971692 161.89283025 181.63594358 201.3790569 221.12217023
# 240.86528356 260.60839689 280.35151021 300.09462354 319.83773687]
graph(
x=x, y=arr, fit=fit,
title='linear (quad)', save_name='linear_quad_curve_fit'
)
2次関数っぽいデータを1次関数でfitすると以下のグラフ。
また、2次関数っぽいデータを2次関数でfitすると以下のコード・グラフ。
def quad(x, a, b, c):
y = a * x ** 2 + b * x + c
return y
par, cov = scipy.optimize.curve_fit(f=quad, xdata=x, ydata=arr)
fit = quad(x, *par)
print(par)
# [ 0.96144676 1.47562492 -0.47895115]
print(cov)
# [[ 2.91900060e-03 -5.54610317e-02 1.66383223e-01]
# [-5.54610317e-02 1.13082139e+00 -3.89336536e+00]
# [ 1.66383223e-01 -3.89336536e+00 1.90009285e+01]]
print(fit)
# [ 0.51600726 1.5202958 2.52315604 3.52458798 4.52459161 5.52316695
# 6.52031398 7.51603271 8.51032314 9.50318526 10.49461909 11.48462461
# 12.47320183 13.46035075 14.44607137 15.43036368 16.41322769 17.3946634
# 18.37467081 19.35324992]
graph(x=x, y=arr, fit=fit, title='quad', save_name='quad_curve_fit')
scipy.optimize.curve_fitでパラメータの誤差を求める
polyfitでパラメータの誤差を求められたように、curve_fitでもパラメータの誤差を求めることが可能。polyfitと同じように共分散行列の対角線の平方根を計算すればいい。
1次関数
# 1次関数 x = np.arange(20) arr = get_poly(power=1) par, cov = scipy.optimize.curve_fit(f=linear, xdata=x, ydata=arr) fit = linear(x, *par) # 対角行列の平方根がパラメータのエラー perr = np.sqrt(np.diag(cov)) print(perr) # [0.01212982 0.13479937]
covは既に取得しているのでそのまま対角線を抽出、平方根をつける。
2次関数
arr = get_poly(power=2) par, cov = scipy.optimize.curve_fit(f=linear, xdata=x, ydata=arr) fit = linear(x, *par) perr = np.sqrt(np.diag(cov)) print(perr) # [ 1.19521437 13.2824833 ]
2次関数も同様。
sin, log, exp関数をfit（scipy.optimize.curve_fit）
上で見てきた多項式函数については既にpolyfitで作成したので驚きはない。ここからはscipy.optimize.curve_fitで多項式以外の関数をfitしてみる。簡単な関数のみを使用するが、やり方は同じなので適宜、関数などを変更して使用していただければと思う。
sin関数っぽいデータをsin関数でfit
def sin_func(x, a, b):
y = a * np.sin(x) + b
return y
def get_sin(x):
# sin関数っぽいデータの作成
arr = []
for num, i in enumerate(x):
np.random.seed(seed=num) # 乱数1
coefficient1 = np.random.rand()
np.random.seed(seed=num + 1) # 乱数2
coefficient2 = np.random.rand() / 10
# 各xで乱数を変えながらsin関数を作成
val = coefficient1 * np.sin(i) + coefficient2
arr.append(val)
arr = np.array(arr)
return arr
# 横軸の値を変更（0から9.98まで0.02刻み）
x = np.arange(0, 10, 0.02)
arr = get_sin(x=x)
まず初めはsin関数だ。fit対象の関数は以下のように2つの変数を持ち、fitするデータはsinぽくなるように乱数を使用してデータを作成した。
$y\ =\ A\ *\ \sin(x)\ +\ B$
fit・グラフ化の方法は多項式と同じようにすればいい。変数が2種類なのでparで返される要素数も2となる。
par, cov = scipy.optimize.curve_fit(f=sin_func, xdata=x, ydata=arr)
fit = sin_func(x, *par)
print(f"par: {par}")
# par: [0.4563479 0.05557946]
graph(x=x, y=arr, fit=fit, title='sin', save_name='sin_curve_fit')
log関数っぽいデータをlog関数でfit
def log_func(x, a):
y = np.log(x) + a
return y
def get_log(x):
# log関数っぽいデータの作成
arr = []
for num, i in enumerate(x):
np.random.seed(seed=num) # 乱数
coefficient1 = np.random.rand() / 10
# 各xで乱数を変えながらlog関数を作成
val = np.log(i) + coefficient1
arr.append(val)
arr = np.array(arr)
return arr
# 横軸の値を変更（1から5 + 0.01まで0.01刻み）
x = np.arange(1, 5 + 0.01, 0.01)
arr = get_log(x=x)
logも同じようにデータを作成してfitする。ただし、logの場合はxが0だと負の無限になるのでスタートを1としておいた。
par, cov = scipy.optimize.curve_fit(f=log_func, xdata=x, ydata=arr)
fit = log_func(x, *par)
print(f"par: {par}")
# par: [0.04786226]
graph(x=x, y=arr, fit=fit, title='log', save_name='log_curve_fit')
exp関数っぽいデータをexp関数でfit
def exp_func(x, a, b):
y = a * np.exp(x) + b
return y
def get_exp(x):
# exp関数っぽいデータの作成
arr = []
for num, i in enumerate(x):
np.random.seed(seed=num)
coefficient1 = np.random.rand()
# 各xで乱数を変えながらexp関数を作成
val = np.exp(i) + coefficient1
arr.append(val)
arr = np.array(arr)
return arr
# 横軸の値を変更（0から1.5まで0.01刻み）
x = np.arange(0, 1.5, 0.001)
arr = get_exp(x=x)
exp関数の場合はxが大きすぎると値がぶっ飛ぶのでセーブするためxを小さめにした。
par, cov = scipy.optimize.curve_fit(f=exp_func, xdata=x, ydata=arr)
fit = exp_func(x, *par)
print(f"par: {par}")
# par: [0.49307665 0.47458626]
graph(x=x, y=arr, fit=fit, title='exp', save_name='exp_curve_fit')
ガウシアンと直線に乗ったガウシアンをfit（scipy.optimize.curve_fit）
お次は、天文系の研究室出身の執筆者がよく使っていたガウシアンをfitしてみる。さらに通常のガウシアンだけではなく、直線に乗ったガウシアンもfitしてみる。
ガウシアンとはつりがね型の関数で、面積を1とした場合は以下の形の式となる。$\sigma,\ \mu$はそれぞれ標準偏差と平均値を表す。
$y\ =\ \cfrac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\cfrac{(x\ -\ \mu)^2}{2\sigma^2} \right\}$
これを直線に乗せるとイメージは以下の式。
$y\ =\ \cfrac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\cfrac{(x\ -\ \mu)^2}{2\sigma^2} \right\}\ +\ ax\ +\ b$
これで斜めになったグラフとなる。グラフ形状は以下で示すので見ていってほしい。
ガウシアンっぽいデータをガウシアンでfit
def gaussian_func(x, sigma, mu, b):
eq1 = 1 / (np.sqrt(2 * np.pi) * sigma)
eq2 = np.exp(-((x - mu) ** 2 / (2 * sigma ** 2)))
y = eq1 * eq2 + b
return y
def get_gaussian(x):
# ガウシアンっぽいデータの作成
arr = []
for num, i in enumerate(x):
np.random.seed(seed=num)
coefficient1 = np.random.rand() / 100
# 各xで乱数を変えながらガウシアン関数を作成
eq1 = 1 / (np.sqrt(2 * np.pi))
eq2 = np.exp(-((i - 1) ** 2) / 2)
val = eq1 * eq2 + coefficient1
arr.append(val)
arr = np.array(arr)
return arr
# 横軸の値を変更（-5から4.995まで0.01刻み）
x = np.arange(-5, 5, 0.005)
arr = get_gaussian(x=x)
ガウシアンと言ってもなんら怖くない。先ほど同様、ガウシアンの関数をdefしてそれをfitするのみ。もし、fitがうまくいかない時はcurve_fitの引数p0でパラメータの初期値を決めることができる。
今回の例では初期値がなくてもいけると思うが、もしうまくfitができない時はp0を設定してみるといいだろう。
# 初期値の設定は引数p0で行う
p0 = (0.1, 0.01, 0.1)
par, cov = scipy.optimize.curve_fit(f=gaussian_func, xdata=x, ydata=arr, p0=p0)
fit = gaussian_func(x, *par)
print(f"par: {par}")
# par: [0.49307665 0.47458626]
graph(x=x, y=arr, fit=fit, title='gaussian', save_name='gaussian_curve_fit')
直線に乗ったガウシアンっぽいデータをfit
def original_func(x, sigma, mu, a, b):
eq1 = 1 / (np.sqrt(2 * np.pi) * sigma)
eq2 = np.exp(-((x - mu) ** 2 / (2 * sigma ** 2)))
y = eq1 * eq2 + a * x + b
return y
def get_original(x):
# 直線に乗ったガウシアンっぽいデータの作成
arr = []
for num, i in enumerate(x):
np.random.seed(seed=num + 1)
coefficient1 = np.random.rand() / 10
np.random.seed(seed=num + 2)
coefficient2 = np.random.rand() / 10
# # 各xで乱数を変えながら直線に乗ったガウシアンを作成
# val = (2 * np.sin(i)) + (coefficient1 * i + coefficient2)
eq1 = 1 / (np.sqrt(2 * np.pi))
eq2 = np.exp(-((i - 1) ** 2) / 2)
val = eq1 * eq2 + coefficient1 * i + coefficient2
arr.append(val)
arr = np.array(arr)
return arr
# 横軸の値を変更（-5から4.995まで0.01刻み）
x = np.arange(-5, 5, 0.005)
arr = get_original(x=x)
直線に乗っていようがいまいが、式として表せるなら簡単にグラフ化できる。
par, cov = scipy.optimize.curve_fit(f=original_func, xdata=x, ydata=arr, )
fit = original_func(x, *par)
print(f"par: {par}")
# par: par: [0.99708121 1.00936314 0.04816927 0.04833364]
graph(x=x, y=arr, fit=fit, title='original', save_name='original_curve_fit')
既存のデータにfit（scipy.optimize.curve_fit）
直線に乗ったガウシアンも計算できた。これで無敵かと思ったけど、実はそうではなかった。大学院終了後も研究を続けているんだけど、その中で既存のデータ配列に他のデータ配列をfitしたい場面が出てきた。
予め関数形がわかっているものではないのでこれは大変。ということで試行錯誤して結局たどり着いたのが、逆算してfitさせるというもの。以下では既存データにfitする方法について解説。
fit対象のデータをarr_target、fitしたいデータをarrとすると以下のイメージでfitすることにする。
$arr\_target\ =\ c * arr\ +\ (ax+\ b)$
この式のarrを元のデータとしてガウシアン、arr_targetを直線に乗ったガウシアンとして進める。
fit対象のデータを作成
def get_gaussian_linear(x):
# 直線にガウシアンが乗ってるっぽいデータの作成
arr = []
for num, i in enumerate(x):
np.random.seed(seed=num)
coefficient1 = np.random.rand() / 10 + 0.5
# 各xで乱数を変えながらガウシアン関数を作成
eq1 = 1 / (np.sqrt(2 * np.pi))
eq2 = np.exp(-((i - 1) ** 2) / 2)
eq3 = 0.05 * i + coefficient1
val = eq1 * eq2 + eq3
arr.append(val)
arr = np.array(arr)
return arr
# 横軸の値を変更（-5から4.995まで0.01刻み）
x = np.arange(-5, 5, 0.005)
# fitの相手プロット
arr_target = get_gaussian_linear(x=x)
# これからfitしたいプロット
arr = get_gaussian(x=x)
まずはfit対象の既存データとfitしたいデータを作成。一応fit対象データの中身は変えてある。get_gaussian_linearでfit対象データを作成し、その配列をarr_targetとして使用。一方でfitしたいデータはarrとした。
fit用の関数と元に戻す関数
def gau_lin_target_func(x, c, a, b):
"""通常のガウシアンを、直線に乗ったガウシアンにfitするための逆算関数
fitする際にのみ使用
Parameters
----------
x : array
横軸の値
c : float
係数1
a : float
係数2
b : float
係数3
Returns
-------
numpy.ndarray
fitした後の関数
"""
numerator = arr_target - (a * x + b) # 分母
denominator = c # 分子
y = numerator / denominator
return y
def gau_lin_func(x, c, a, b):
"""fitパラメータを引数c, a, bに入れることで
通常のガウシアンを直線に乗ったガウシアンにfitする
Parameters
----------
x : array
横軸の値
c : float
係数1
a : float
係数2
b : float
係数3
Returns
-------
numpy.ndarray
通常のガウシアンを直線に乗ったガウシアンにfitした後の配列
"""
y = c * arr + a * x + b
return y
ここが今回の章のキモ。fitしたい数式を逆算した関数gau_lin_target_funcでデータをfitする。こうすることで、既存の知られた関数じゃなくてもfitすることができる。
fitしたあとはパラメータを得ることができるので、このパラメータを本来の関数gau_lin_funcに代入することでfitが完了する。
これを見つけた時は感動したよ。で、この記事を書こうとした。fitは今まで通り。
par, cov = scipy.optimize.curve_fit(f=gau_lin_target_func, xdata=x, ydata=arr)
fit = gau_lin_func(x, *par)
print(f"par: {par}")
# par: [1.04663899 0.04940854 0.03884187]
グラフ化は元のデータ、fit対象、fit後の元のデータの3種類なので新規で関数を作成した。