こんな人にオススメ
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のグラフ即席作成コード
続きを見る
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種類なので新規で関数を作成した。
fitでデータの傾向がわかる
今回はpythonを使用してグラフをfitした。fit対象はよく知られた関数、既存のデータの2種類でfitした。さらにパラメータの誤差も計算した。
実際にデータをfitすることでそのデータの傾向がよりわかりやすくなる。2次関数だと思ってfitしたら全然形状が異なることもある。実際、研究で多項式だと思ったデータが全然fitできなかったこともある。
既存のデータにfitする方法は色々と調べても出なかったので、本記事で方法を知って活用していただければ幸いだ。
関連記事
-
-
【python&関数化】defとかargsとかを使って関数を作成する
続きを見る
-
-
【python3&zip関数】python3系にてzipは一回使ったら消える
続きを見る
-
-
【max(), sorted()&key】max関数とかで使うkeyを活用
続きを見る
-
-
【Python3.10&パターンマッチ】match文を試してみる
続きを見る