カテゴリー

Python基礎

【python&フィッティング】polyfitとcurve_fitでfitting

2021年7月15日

こんな人にオススメ

pythonを使ってデータを関数とか既存の他のデータにフィットさせたい(フィッティング)。どうしたらいい?

ということで今回はpythonを使ってデータをフィッティングする方法について解説する。pythonでのフィッティングの方法はいくつかあるが、今回はnumpypolyfitscipycurvefitにフォーカスする。

python環境は以下。

  • Python 3.9.4
  • numpy 1.20.3
  • scipy 1.6.3
  • plotly 4.14.3
  • plotly-orca 3.4.2
スポンサーリンク
スポンサーリンク

運営者のメガネと申します。TwitterInstagramも運営中。

自己紹介はこちらから、お問い合わせはこちらからお願いいたします。

運営者メガネ

下準備

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

まずは今回使用するモジュールのimportplotly_layout_templateは自作関数で、2つ上のディレクトリに置いているので../../で移動している。詳細は以下参照。

【随時更新 備忘録】plotlyのグラフ即席作成コード

続きを見る

1, 2次関数をfit(np.polyfit

まずはnumpypolyfit。これは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_tff_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の次数を表すorder2にするというとこ。これに伴ってパラメータは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_tfTrueにして共分散行列を出力。そのあと、その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したい関数と横・縦軸の値を指定すれば良い。戻り値はパラメータと共分散行列。それぞれparcovに代入した。

あとは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の場合はx0だと負の無限になるのでスタートを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文を試してみる

続きを見る

関連コンテンツ

スポンサーリンク

Amazonのお買い物で損したない人へ

1回のチャージ金額通常会員プライム会員
¥90,000〜2.0%2.5%
¥40,000〜1.5%2.0%
¥20,000〜1.0%1.5%
¥5,000〜0.5%1.0%

Amazonギフト券にチャージすることでお得にお買い物できる。通常のAmazon会員なら最大2.0%、プライム会員なら2.5%還元なのでバカにならない。

ゲットしたポイントは通常のAmazonでのお買い物に使えるからお得だ。一度チャージしてしまえば、好きなタイミングでお買いものできる。

なお、有効期限は10年だから安心だ。いつでも気軽にAmazonでお買い物できる。

Amazonチャージはここから出来るで

もっとお得なAmazon Prime会員はこちらから

30日間無料登録

執筆者も便利に使わせてもらってる

スポンサーリンク

  • この記事を書いた人

メガネ

独学でpythonを学び天文学系の大学院を修了。 ガジェット好きでMac×Android使い。色んなスマホやイヤホンを購入したいけどお金がなさすぎて困窮中。 元々、人見知りで根暗だったけど、人生楽しもうと思って良い方向に狂ったために今も人生めちゃくちゃ楽しい。 pythonとガジェットをメインにブログを書いていますので、興味を持たれましたらちょこちょこ訪問してくだされば幸いです🥰。 自己紹介→変わって楽しいの繰り返し

-Python基礎
-, , ,