Pythonでフィッティング (線形/非線形/多変数を例示し、誤差評価まで)

2020 7/05

2020.07.05 2020.07.05

今回はpythonを使ったフィッティング(fitting)のやり方をご紹介します。

実験や観測で、直線っぽい、二次曲線っぽい、2変数のガウシアンっぽいデータが得られた時、機械学習なんて使わずに最小二乗法を使ったfittingを行うと大体データを数式で表せます。

pythonでfittingを行うのは、scipyのcurve_fitを使えば、かなり簡単です。

得られた数式の誤差評価(今回紹介するのは、カイ二乗評価と決定係数評価)もpythonなら数行で計算してくれます。

下記では、最小二乗法を使用したfittingの原理や評価方法の概念の紹介を行い、線形fittingを例にpythonでの実装方法を示した後に、多変数の場合も2Dガウシアンを例に説明していきます。同時にfittingの良し悪しを判断できるカイ二乗評価と決定係数評価も実装していきます。

「Fitting」とは

Fittingは、日本語にすると「曲線あてはめ」などとなります。

測定や観測したデータに対して、モデル式のパラメータを変えていき、最もデータに合う(fitする)パラメータ値を決定するのがfittingです。

Fittingの何が嬉しいのかというと、得られたモデル式を使って、実測データがない部分の値を推定(内挿、外挿)できるのです。例えば1日の気温を測定した時に、7時(20℃)と9時(22℃)のデータがあれば、8時の気温が21℃だと推定できます。

fittingの場合、機械学習とは異なり、モデル式は自分で考える必要があります。

モデル式は、理論的に「この式だ！」という式がある場合もあれば、よくわからないけど直線ぽいからとりあえず直線の式で当てはめてみるとか、全くどんな式か見当もつかないからべき乗関数式(テーラー展開)でやってみる、など様々です。

イメージとしては次の図のようになります。

Fitting曲線(モデル式)内のパラメータ$ a, b, c, d $を変化させて、実測データとモデル式の差が最も少なくなるようなパラメータセットを決定するのが最小二乗法によるfittingです。

ここで「実測データとモデル式の差」とは、次の図中の赤い線の長さの合計(残差平方和)のことです。

この差が少なければ確かに、モデル式と実測データは合っていて(fitしていて)、モデル式を使えば実測データがない部分での$ y $の予測もできそうなのがわかります。

Fittingの良し悪しの評価方法

Fitting後に、どれだけモデル式がデータに合ったかを定量的に示す方法を2つ紹介します。

カイ二乗

データ値$ y_o $と、モデル式から計算される予測値$ y_e $と、を使って次の「カイ二乗」を計算します。

$$ \chi^2 = \sum\frac{(y_o-y_e)^2}{y_e} $$

この数値は、「二乗」と言いつつ負も許される値ですが、これが0に近ければ近いほどに、モデル式がデータ値に合っている(fitしている)ことになります。

色々と難しい数式を追うと、この$ \chi^2 $の値は、「自由度」の数$ d $と同程度か小さければ、モデル式の精度がよく、「自由度」の数$ d $よりも非常に大きければ精度が悪いです。

$$ \chi^2 \leq d \Rightarrow モデル式の精度高い \\ \chi^2 \gg d \Rightarrow モデル式の精度低い $$

「自由度」とは、「データ数 – モデル式のパラメータ数」です。前節の図の例だと、データの数は16で、パラメータの数は$ a, b, c, d $の4つなので、自由度は12です。

カイ二乗評価の枠組み内には実はもう一つ「p値(p-value)」というものがあり、次のような式で表されます。

$$ p = \frac{2}{2^{d/2}\Gamma\left(\frac{d}{2}\right)}\int_{\chi_0}^{\infty}t^{d-1}e^{-\frac{t^2}{2}}{\rm d}t $$

ここで$ \chi_0 $は「カイ二乗」の平方根で、$ \Gamma(z) $はガンマ関数で、$ d $は自由度です。

このp値は0(=0%)から1(=100%)の間の数値で、モデル式がデータにfitする確率のようなものです。

$$ p\simeq 1 \Rightarrow モデル式の精度高い \\ p\simeq 0 \Rightarrow モデル式の精度低い $$

ここまでガンマ関数が出てきたりと、かなり数学が濃くなってしまいましたが、pythonでこれを計算する時には、何も考えずに次のようにコードを書き込めばい良いだけです。pythonバンザイ。

#カイ二乗計算例
from scipy import stats

o = [1,2,3,4,5,6] #実測データ例
e = [2,3,3,4,5,6] #モデル計算された数値例
chi2 = stats.chisquare(o, f_exp = e)

print(chi2)
# -> Power_divergenceResult(statistic=0.8333333333333333, pvalue=0.9748586976330175)

#カイ二乗計算例

from scipy import stats

o = [1,2,3,4,5,6] #実測データ例

e = [2,3,3,4,5,6] #モデル計算された数値例

chi2 = stats.chisquare(o, f_exp = e)

print(chi2)

# -> Power_divergenceResult(statistic=0.8333333333333333, pvalue=0.9748586976330175)

ここでは自由度を入力していませんが、その場合には「データ数 – 1」が自由度として設定されます。

結果は変数「chi2」に格納されて、chi2[0]がカイ二乗値、chi2[1]がp値となります。

自由度設定などを詳しく知りたい方はここをご参照ください。

決定係数

ここではもう一つの「fittingの良し悪しを評価できる数値」として、決定係数( $ R^2 $ )を紹介します。

決定係数は次のような式で定義されています。

$$ R^2 = 1- \frac{\sum \left( y_o – y_e \right)^2}{\sum \left( y_o – \overline{y_o} \right)^2} $$

ここで$ y_o $ は測定値で、$ y_e $ はモデル式による予測値です。$ \overline{y_o} $ は観測値の平均を表しています。

この数値はイメージがしやすいですね。

$$ R^2 \simeq 1 \Rightarrow モデル式の精度高い \\ R^2 \simeq 0 \Rightarrow モデル式の精度低い $$

実は、この$ R^2 $ も二乗がついているものの、負の値を取り得ます。もちろん1から離れれば離れるほどに、モデル式の精度が悪いことを表しています。

Pythonでこれを計算する時には、単純に次のように計算させます。

#決定係数の計算
import numpy as np

o = np.array([1,2,3,4,5,6]) #実測データ例
e = np.array([2,3,3,4,5,6]) #モデル計算された数値例

residuals =  o - e #残渣
rss = np.sum(residuals**2)      #残差平方和: residual sum of squares = rss
tss = np.sum((o-np.mean(o))**2) #全平方和: total sum of squares = tss
r_squared = 1 - (rss / tss)     #決定係数R^2

print(r_squared)
# -> 0.8857142857142857

#決定係数の計算

import numpy as np

o = np.array([1,2,3,4,5,6]) #実測データ例

e = np.array([2,3,3,4,5,6]) #モデル計算された数値例

residuals = o - e #残渣

rss = np.sum(residuals**2) #残差平方和: residual sum of squares = rss

tss = np.sum((o-np.mean(o))**2) #全平方和: total sum of squares = tss

r_squared = 1 - (rss / tss) #決定係数R^2

print(r_squared)

# -> 0.8857142857142857

Pythonでfittingする方法 <線形関数を例に>

これはかなり簡単で、scipyのcurve_fitというのを使用します。

メインの計算は次の1行でできます。

#scipyのcurve_fitを使用したfittingのメインの部分
popt, pcov = curve_fit(linear_func,x_observed,y_observed) #poptは最適推定値、pcovは共分散が出力される
perr = np.sqrt(np.diag(pcov)) #推定されたパラメータの各々の誤差

#scipyのcurve_fitを使用したfittingのメインの部分

popt, pcov = curve_fit(linear_func,x_observed,y_observed) #poptは最適推定値、pcovは共分散が出力される

perr = np.sqrt(np.diag(pcov)) #推定されたパラメータの各々の誤差

linear_funcは定義された関数(モデル式)で、x_observedとy_observedが観測データです。

Fittingで得られたパラメータはpoptに格納されますが、その推定誤差の情報はpcovに格納されます。

pcovは実際には共分散行列なので、上記コードのように対角成分の平方根をとることで、推定誤差を算出します。

実際に線形関数をモデル式にした場合の例を示します。

#linear_fitting_0.py

import numpy as np
from scipy.optimize import curve_fit

#fittingを行うメインの関数
def fitting_linear(data):
    #モデル式を表す関数
    def linear_func(X, a, b): # １次式
        Y = a * X + b
        return Y

    x_observed = data["x"]
    y_observed = data["y"]

    #fittingのメイン計算部分
    popt, pcov = curve_fit(linear_func,x_observed,y_observed) #poptは最適推定値、pcovは共分散が出力される
    perr = np.sqrt(np.diag(pcov)) #推定されたパラメータの各々の誤差

    #fittingの結果をターミナルに表示
    print("*Result************")
    print("y = ax + b")
    print("a = ", popt[0], "+-", perr[0])
    print("b = ", popt[1], "+-", perr[1])
    print("*******************")

    return 0

if __name__=="__main__":

    #入力データ(観測データ)
    data={  "x": np.array([0, 15, 30, 45, 60, 75, 90, 115, 130, 145, 160, 175, 190, 215, 230, 255, 275, 290, 315, 330, 340, 360, 370]),
            "y": np.array([170, 270, 167, 200, 160, 130, 180, 190, 120, 70, 80, 50, 80, 34, 120, 190, 73, 80, 30, 50, 20, 40, 90])
            }

    #fittingを行う関数の読み出し
    fitting_linear(data)

# ->
# *Result************
# y = ax + b
# a =  -0.4221020913057208 +- 0.08333822145540422
# b =  190.0456436906002 +- 18.04262657531871
# *******************

#linear_fitting_0.py

import numpy as np

from scipy.optimize import curve_fit

#fittingを行うメインの関数

def fitting_linear(data):

#モデル式を表す関数

def linear_func(X, a, b): # １次式

Y = a * X + b

return Y

x_observed = data["x"]

y_observed = data["y"]

#fittingのメイン計算部分

popt, pcov = curve_fit(linear_func,x_observed,y_observed) #poptは最適推定値、pcovは共分散が出力される

perr = np.sqrt(np.diag(pcov)) #推定されたパラメータの各々の誤差

#fittingの結果をターミナルに表示

print("*Result************")

print("y = ax + b")

print("a = ", popt[0], "+-", perr[0])

print("b = ", popt[1], "+-", perr[1])

print("*******************")

return 0

if __name__=="__main__":

#入力データ(観測データ)

data={ "x": np.array([0, 15, 30, 45, 60, 75, 90, 115, 130, 145, 160, 175, 190, 215, 230, 255, 275, 290, 315, 330, 340, 360, 370]),

"y": np.array([170, 270, 167, 200, 160, 130, 180, 190, 120, 70, 80, 50, 80, 34, 120, 190, 73, 80, 30, 50, 20, 40, 90])

}

#fittingを行う関数の読み出し

fitting_linear(data)

# ->

# *Result************

# y = ax + b

# a = -0.4221020913057208 +- 0.08333822145540422

# b = 190.0456436906002 +- 18.04262657531871

# *******************

このコードの書き方では、関数をまず定義した後に、if name=="main":の部分から処理が開始されます。

観測データがnumpy配列として「data」に格納されて、fitting_linearに渡されます。

fitting_linear関数内では、モデル式として、linear_funcが定義されていて、curve_fitを使ってfittingが行われています。超簡単！

curve_fitでは、最初の入力変数がモデル式で、2番目がx(独立変数、説明変数)で、3番目がy(従属変数、被説明変数)となっています。

(このxは変数の数としては1つしか入力できないのですが、中身は配列でもよいので、多変数にできます。詳細は後ほど。)

評価値計算とグラフ表示も含めた例示

先述のカイ二乗、p値、決定係数の計算と、測定データおよびfitting結果のグラフ化まで含めたコードを示しておきます。

#linear_fitting.py

import numpy as np
from scipy.optimize import curve_fit
from scipy import stats
from matplotlib import pyplot as plt

#fittingを行うメインの関数
def fitting_linear(data):
    #モデル式を表す関数
    def linear_func(X, a, b): # １次式
        Y = a * X + b
        return Y

    #fitting結果のグラフ化関数
    def plot_fit_result(data, fit_result):
        plt.scatter(data["x"], data["y"], c="blue", s=20)   #実測データ値は散布図でplot
        plt.plot(fit_result["x"], fit_result["y"], c="red") #fitting結果は線でplot
        plt.show(block=False)                               #グラフを表示
        plt.savefig("linear_fitting_result.png", dpi=300)    #グラフを画像データで保存
        return 0

    x_observed = data["x"]
    y_observed = data["y"]

    #fittingのメイン計算部分
    popt, pcov = curve_fit(linear_func,x_observed,y_observed) #poptは最適推定値、pcovは共分散が出力される
    perr = np.sqrt(np.diag(pcov)) #推定されたパラメータの各々の誤差

    #観測データと推定データの間のカイ自乗を計算
    o = y_observed #観測データ
    e = linear_func(x_observed, popt[0], popt[1]) #推定データ
    chi2 = stats.chisquare(o, f_exp = e) #カイ自乗計算のメイン部分。chi2には[カイ二乗, p値]の2つが出力される。

    #決定係数(R二乗)計算
    residuals =  o - e #残渣
    rss = np.sum(residuals**2)      #残差平方和: residual sum of squares = rss
    tss = np.sum((o-np.mean(o))**2) #全平方和: total sum of squares = tss
    r_squared = 1 - (rss / tss)

    #fittingの結果をターミナルに表示
    print("*Result************")
    print("y = ax + b")
    print("a = ", popt[0], "+-", perr[0])
    print("b = ", popt[1], "+-", perr[1])
    print("X-squared = ", chi2[0])
    print("p-value = ", chi2[1])
    print("R^2 = ", r_squared)
    print("*******************")
    statistics_numbers = {  "X-squared": chi2[0],
                            "p-value": chi2[1],
                            "R^2": r_squared}

    #グラフ用に、fittingの結果を示す線を作成
    fit_x = np.linspace(min(data["x"]), max(data["x"]), 100)
    fit_y = linear_func(fit_x, popt[0], popt[1])
    fit_result={"x":fit_x, "y":fit_y}

    #グラフでfitting結果を表示
    plot_fit_result(data, fit_result)

    return popt[0], popt[1], statistics_numbers #fittingで得られたパラメータと統計数値を返す。

if __name__=="__main__":

    #入力データ(観測データ)
    data={  "x": np.array([0, 15, 30, 45, 60, 75, 90, 115, 130, 145, 160, 175, 190, 215, 230, 255, 275, 290, 315, 330, 340, 360, 370]),
            "y": np.array([170, 270, 167, 200, 160, 130, 180, 190, 120, 70, 80, 50, 80, 34, 120, 190, 73, 80, 30, 50, 20, 40, 90])
            }

    #fittingを行う関数の読み出し
    a, b, statistics_numbers = fitting_linear(data)

    #グラフ表示のための入力待機
    Answer="n"
    while Answer != "y":
        print("Do you want to stop this program? (y/n)")
        Answer = input(">> ")
    print("End")

# *Result************
# y = ax + b
# a =  -0.4221020913057208 +- 0.08333822145540422
# b =  190.0456436906002 +- 18.04262657531871
# X-squared =  478.9130766425304
# p-value =  1.8042438803936564e-87
# R^2 =  0.5498729579176961
# *******************

#linear_fitting.py

import numpy as np

from scipy.optimize import curve_fit

from scipy import stats

from matplotlib import pyplot as plt

#fittingを行うメインの関数

def fitting_linear(data):

#モデル式を表す関数

def linear_func(X, a, b): # １次式

Y = a * X + b

return Y

#fitting結果のグラフ化関数

def plot_fit_result(data, fit_result):

plt.scatter(data["x"], data["y"], c="blue", s=20) #実測データ値は散布図でplot

plt.plot(fit_result["x"], fit_result["y"], c="red") #fitting結果は線でplot

plt.show(block=False) #グラフを表示

plt.savefig("linear_fitting_result.png", dpi=300) #グラフを画像データで保存

return 0

x_observed = data["x"]

y_observed = data["y"]

#fittingのメイン計算部分

popt, pcov = curve_fit(linear_func,x_observed,y_observed) #poptは最適推定値、pcovは共分散が出力される

perr = np.sqrt(np.diag(pcov)) #推定されたパラメータの各々の誤差

#観測データと推定データの間のカイ自乗を計算

o = y_observed #観測データ

e = linear_func(x_observed, popt[0], popt[1]) #推定データ

chi2 = stats.chisquare(o, f_exp = e) #カイ自乗計算のメイン部分。chi2には[カイ二乗, p値]の2つが出力される。

#決定係数(R二乗)計算

residuals = o - e #残渣

rss = np.sum(residuals**2) #残差平方和: residual sum of squares = rss

tss = np.sum((o-np.mean(o))**2) #全平方和: total sum of squares = tss

r_squared = 1 - (rss / tss)

#fittingの結果をターミナルに表示

print("*Result************")

print("y = ax + b")

print("a = ", popt[0], "+-", perr[0])

print("b = ", popt[1], "+-", perr[1])

print("X-squared = ", chi2[0])

print("p-value = ", chi2[1])

print("R^2 = ", r_squared)

print("*******************")

statistics_numbers = { "X-squared": chi2[0],

"p-value": chi2[1],

"R^2": r_squared}

#グラフ用に、fittingの結果を示す線を作成

fit_x = np.linspace(min(data["x"]), max(data["x"]), 100)

fit_y = linear_func(fit_x, popt[0], popt[1])

fit_result={"x":fit_x, "y":fit_y}

#グラフでfitting結果を表示

plot_fit_result(data, fit_result)

return popt[0], popt[1], statistics_numbers #fittingで得られたパラメータと統計数値を返す。

if __name__=="__main__":

#入力データ(観測データ)

data={ "x": np.array([0, 15, 30, 45, 60, 75, 90, 115, 130, 145, 160, 175, 190, 215, 230, 255, 275, 290, 315, 330, 340, 360, 370]),

"y": np.array([170, 270, 167, 200, 160, 130, 180, 190, 120, 70, 80, 50, 80, 34, 120, 190, 73, 80, 30, 50, 20, 40, 90])

}

#fittingを行う関数の読み出し

a, b, statistics_numbers = fitting_linear(data)

#グラフ表示のための入力待機

Answer="n"

while Answer != "y":

print("Do you want to stop this program? (y/n)")

Answer = input(">> ")

print("End")

# *Result************

# y = ax + b

# a = -0.4221020913057208 +- 0.08333822145540422

# b = 190.0456436906002 +- 18.04262657531871

# X-squared = 478.9130766425304

# p-value = 1.8042438803936564e-87

# R^2 = 0.5498729579176961

# *******************

結果のグラフは次のようになります。

モデル式(linear_func部分)を適当に変えれば、非線形でもなんでもモデル式を作成できます。

多変数のfittingを行う場合

上記の例ではyはxのみの関数でした。つまり独立変数(説明変数)は1つでした。

ここでは独立変数が複数の時にはどのようにfittingをすれば良いかを、2変数を例に説明していきます。

ちなみに、モデル式は2次元のガウシアンです。

$$ z = A\frac{1}{\sqrt{2\pi \sigma_x^2}} e^{-\frac{(x-\mu_x)^2}{2\sigma_x^2}} \frac{1}{\sqrt{2\pi \sigma_y^2}} e^{-\frac{(y-\mu_y)^2}{2\sigma_y^2}} $$

1変数と特に変わったところは、curve_fitへの入力とモデル式を表す関数の入力の仕方です。

curve_fitへの入力は次のようにします。

#scipyのcurve_fitを使用したfittingのメインの部分
#2変数の場合
popt, pcov = curve_fit(linear_func, (x_observed, y_observed), z_observed) #poptは最適推定値、pcovは共分散が出力される

#scipyのcurve_fitを使用したfittingのメインの部分

#2変数の場合

popt, pcov = curve_fit(linear_func, (x_observed, y_observed), z_observed) #poptは最適推定値、pcovは共分散が出力される

独立変数(説明変数)は2つをカッコにくくり、1つの変数として入力します。

モデル式の定義では次のようにします。

#2変数のfittingにおけるモデル式定義関数
def two_D_gauss(X, A, sigma_x, sigma_y, mu_x, mu_y): # 2D Gaussian
    x, y = X #説明変数(独立変数)を分離する。
    z = A * 1/(np.sqrt(2*np.pi*sigma_x**2))*np.exp(-(x-mu_x)**2/(2*sigma_x**2)) * 1/(np.sqrt(2*np.pi*sigma_y**2))*np.exp(-(y-mu_y)**2/(2*sigma_y**2))

    return z

#2変数のfittingにおけるモデル式定義関数

def two_D_gauss(X, A, sigma_x, sigma_y, mu_x, mu_y): # 2D Gaussian

x, y = X #説明変数(独立変数)を分離する。

z = A * 1/(np.sqrt(2*np.pi*sigma_x**2))*np.exp(-(x-mu_x)**2/(2*sigma_x**2)) * 1/(np.sqrt(2*np.pi*sigma_y**2))*np.exp(-(y-mu_y)**2/(2*sigma_y**2))

return z

このようにすることで、入力は1変数Xにして、関数中で2つ(=x, y)に分けてあげます。

多変数の場合のfittingコード例 <評価値計算とグラフ表示を含む>

コードは次のようになります。

#twoD_Gauss_fitting.py

import numpy as np
from scipy.optimize import curve_fit
from scipy import stats
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def fitting_2DGauss(data):
    def two_D_gauss(X, A, sigma_x, sigma_y, mu_x, mu_y): # 2D Gaussian
        x, y = X #説明変数(独立変数)を分離する。
        z = A * 1/(np.sqrt(2*np.pi*sigma_x**2))*np.exp(-(x-mu_x)**2/(2*sigma_x**2)) * 1/(np.sqrt(2*np.pi*sigma_y**2))*np.exp(-(y-mu_y)**2/(2*sigma_y**2))

        return z

    def plot_fit_result(data, fit_result):
        #グラフの枠を作っていく
        fig = plt.figure()
        ax = Axes3D(fig)

        #描画
        ax.plot(data["x"],data["y"],data["z"], ms=3, marker="o",linestyle='None', c="blue")         #実測データ値は散布図でplot
        ax.plot_wireframe(fit_result["x"],fit_result["y"],fit_result["z"], rstride=10, cstride=10)  #fitting結果は面(ワイヤーフレーム)でplot
        plt.show(block=False)                                                                       #グラフを表示
        plt.savefig("twoD_Gauss_fitting_result.png", dpi=300)                                       #グラフを画像データで保存
        return 0

    x_observed = data["x"]
    y_observed = data["y"]
    z_observed = data["z"]

    #fittingのメイン計算部分
    popt, pcov = curve_fit(two_D_gauss, (x_observed, y_observed), z_observed) #poptは最適推定値、pcovは共分散が出力される
    perr = np.sqrt(np.diag(pcov)) #推定されたパラメータの各々の誤差

    #Chi2 contingency
    o = z_observed #観測データ
    e = two_D_gauss((x_observed, y_observed), popt[0], popt[1], popt[2], popt[3], popt[4]) #推定データ
    chi2 = stats.chisquare(o, f_exp = e) #カイ自乗計算のメイン部分。chi2には[カイ二乗, p値]の2つが出力される。

    #R2 calc
    residuals =  o - e #残渣
    rss = np.sum(residuals**2)      #残差平方和: residual sum of squares = rss
    tss = np.sum((o-np.mean(o))**2) #全平方和: total sum of squares = tss
    r_squared = 1 - (rss / tss)     #決定係数R^2

    #fittingの結果をターミナルに表示
    print("*Result************")
    print("z = A * 1/sqrt(2*pi*sigma_x^2) * exp(-(x-mu_x)^2/(2*sigma_x^2)) * 1/sqrt(2*pi*sigma_y^2) * exp(-(y-mu_y)^2/(2*sigma_y^2))")
    print("A = ", popt[0], "+-", perr[0])
    print("sigma_x = ", popt[1], "+-", perr[1])
    print("sigma_y = ", popt[2], "+-", perr[2])
    print("mu_x = ", popt[3], "+-", perr[3])
    print("mu_y = ", popt[4], "+-", perr[4])
    print("X-squared = ", chi2[0])
    print("p-value = ", chi2[1])
    print("R^2 = ", r_squared)
    print("*******************")
    statistics_numbers = {  "X-squared": chi2[0],
                            "p-value": chi2[1],
                            "R^2": r_squared}

    #グラフ用に、fittingの結果を示す曲面を作成
    fit_x = np.linspace(min(data["x"]), max(data["x"]), 100)
    fit_y = np.linspace(min(data["y"]), max(data["y"]), 100)
    X, Y = np.meshgrid(fit_x, fit_y)
    fit_z = two_D_gauss((X, Y), popt[0], popt[1], popt[2], popt[3], popt[4])
    fit_result={"x":X, "y":Y, "z":fit_z}

    #グラフでfitting結果を表示
    plot_fit_result(data, fit_result)

    return popt[0], popt[1], popt[2], popt[3], popt[4], statistics_numbers

if __name__=="__main__":

    #入力データ(観測データ)
    data={  "x": np.array([8.58889267, 3.72711155, 5.55128778, 9.55656549, 7.36669598, 8.16205139, 1.0108656,  9.2848807,  6.09109169, 5.96553436, 0.91784135, 3.45186243, 6.62752523, 4.41713488, 5.51487786, 7.03712491, 5.89401231, 0.49932759, 5.6179184,  7.66358472, 9.1090833,  0.92909947, 9.0252139,  4.60960406, 4.5201847,  9.99425487, 1.6242374, 7.09370584, 1.60624077, 8.10776768]),
            "y": np.array([0.35147175, 5.34886729, 1.66500123, 3.0841038,  0.45062413, 2.38576126,6.7483453,  7.8238275,  6.95201635, 3.28954449, 4.94031866, 5.24121363,2.98541245, 4.6310814 , 9.84784293, 5.01134921, 3.98072447, 7.27905317,8.6333097 , 0.26169537, 2.90017181, 7.89069194, 4.57119669, 0.0692848,4.19335457, 3.30674764, 6.04152128, 3.24620837, 9.81251081, 5.88231951]),
            "z": np.array([0.02602832, 0.16052722, 0.09373329, 0.04550639, 0.04102718, 0.06941357,0.06164272, 0.0409449 , 0.13394293, 0.14272363, 0.07005369, 0.15429319,0.1218274 , 0.17222531, 0.04722236, 0.1404263 , 0.15967134, 0.04300611,0.08314955, 0.03425328, 0.05417651, 0.04427107, 0.07115732, 0.04542704,0.16839289, 0.03772248, 0.08840031, 0.11684266, 0.0257566,  0.09902988])
            }

    #fittingを行う関数の読み出し
    A, sigama_x, sigma_y, mu_x, mu_y, statistics_numbers = fitting_2DGauss(data)

    #入力待機
    Answer="n"
    while Answer != "y":
        print("Do you want to stop this program? (y/n)")
        Answer = input(">> ")
    print("End")

# *Result************
# z = A * 1/sqrt(2*pi*sigma_x^2) * exp(-(x-mu_x)^2/(2*sigma_x^2)) * 1/sqrt(2*pi*sigma_y^2) * exp(-(y-mu_y)^2/(2*sigma_y^2))
# A =  9.99999979661869 +- 1.0783043023518633e-07
# sigma_x =  2.999999985388892 +- 2.8280927320654662e-08
# sigma_y =  2.999999956090755 +- 3.205619837828646e-08
# mu_x =  4.999999999911432 +- 2.72673909345384e-08
# mu_y =  5.000000016625016 +- 3.18675041452599e-08
# X-squared =  3.187660457301885e-15
# p-value =  1.0
# R^2 =  0.9999999999999969
# *******************

100

101

102

103

#twoD_Gauss_fitting.py

import numpy as np

from scipy.optimize import curve_fit

from scipy import stats

from matplotlib import pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

def fitting_2DGauss(data):

def two_D_gauss(X, A, sigma_x, sigma_y, mu_x, mu_y): # 2D Gaussian

x, y = X #説明変数(独立変数)を分離する。

z = A * 1/(np.sqrt(2*np.pi*sigma_x**2))*np.exp(-(x-mu_x)**2/(2*sigma_x**2)) * 1/(np.sqrt(2*np.pi*sigma_y**2))*np.exp(-(y-mu_y)**2/(2*sigma_y**2))

return z

def plot_fit_result(data, fit_result):

#グラフの枠を作っていく

fig = plt.figure()

ax = Axes3D(fig)

#描画

ax.plot(data["x"],data["y"],data["z"], ms=3, marker="o",linestyle='None', c="blue") #実測データ値は散布図でplot

ax.plot_wireframe(fit_result["x"],fit_result["y"],fit_result["z"], rstride=10, cstride=10) #fitting結果は面(ワイヤーフレーム)でplot

plt.show(block=False) #グラフを表示

plt.savefig("twoD_Gauss_fitting_result.png", dpi=300) #グラフを画像データで保存

return 0

x_observed = data["x"]

y_observed = data["y"]

z_observed = data["z"]

#fittingのメイン計算部分

popt, pcov = curve_fit(two_D_gauss, (x_observed, y_observed), z_observed) #poptは最適推定値、pcovは共分散が出力される

perr = np.sqrt(np.diag(pcov)) #推定されたパラメータの各々の誤差

#Chi2 contingency

o = z_observed #観測データ

e = two_D_gauss((x_observed, y_observed), popt[0], popt[1], popt[2], popt[3], popt[4]) #推定データ

chi2 = stats.chisquare(o, f_exp = e) #カイ自乗計算のメイン部分。chi2には[カイ二乗, p値]の2つが出力される。

#R2 calc

residuals = o - e #残渣

rss = np.sum(residuals**2) #残差平方和: residual sum of squares = rss

tss = np.sum((o-np.mean(o))**2) #全平方和: total sum of squares = tss

r_squared = 1 - (rss / tss) #決定係数R^2

#fittingの結果をターミナルに表示

print("*Result************")

print("z = A * 1/sqrt(2*pi*sigma_x^2) * exp(-(x-mu_x)^2/(2*sigma_x^2)) * 1/sqrt(2*pi*sigma_y^2) * exp(-(y-mu_y)^2/(2*sigma_y^2))")

print("A = ", popt[0], "+-", perr[0])

print("sigma_x = ", popt[1], "+-", perr[1])

print("sigma_y = ", popt[2], "+-", perr[2])

print("mu_x = ", popt[3], "+-", perr[3])

print("mu_y = ", popt[4], "+-", perr[4])

print("X-squared = ", chi2[0])

print("p-value = ", chi2[1])

print("R^2 = ", r_squared)

print("*******************")

statistics_numbers = { "X-squared": chi2[0],

"p-value": chi2[1],

"R^2": r_squared}

#グラフ用に、fittingの結果を示す曲面を作成

fit_x = np.linspace(min(data["x"]), max(data["x"]), 100)

fit_y = np.linspace(min(data["y"]), max(data["y"]), 100)

X, Y = np.meshgrid(fit_x, fit_y)

fit_z = two_D_gauss((X, Y), popt[0], popt[1], popt[2], popt[3], popt[4])

fit_result={"x":X, "y":Y, "z":fit_z}

#グラフでfitting結果を表示

plot_fit_result(data, fit_result)

return popt[0], popt[1], popt[2], popt[3], popt[4], statistics_numbers

if __name__=="__main__":

#入力データ(観測データ)

data={ "x": np.array([8.58889267, 3.72711155, 5.55128778, 9.55656549, 7.36669598, 8.16205139, 1.0108656, 9.2848807, 6.09109169, 5.96553436, 0.91784135, 3.45186243, 6.62752523, 4.41713488, 5.51487786, 7.03712491, 5.89401231, 0.49932759, 5.6179184, 7.66358472, 9.1090833, 0.92909947, 9.0252139, 4.60960406, 4.5201847, 9.99425487, 1.6242374, 7.09370584, 1.60624077, 8.10776768]),

"y": np.array([0.35147175, 5.34886729, 1.66500123, 3.0841038, 0.45062413, 2.38576126,6.7483453, 7.8238275, 6.95201635, 3.28954449, 4.94031866, 5.24121363,2.98541245, 4.6310814 , 9.84784293, 5.01134921, 3.98072447, 7.27905317,8.6333097 , 0.26169537, 2.90017181, 7.89069194, 4.57119669, 0.0692848,4.19335457, 3.30674764, 6.04152128, 3.24620837, 9.81251081, 5.88231951]),

"z": np.array([0.02602832, 0.16052722, 0.09373329, 0.04550639, 0.04102718, 0.06941357,0.06164272, 0.0409449 , 0.13394293, 0.14272363, 0.07005369, 0.15429319,0.1218274 , 0.17222531, 0.04722236, 0.1404263 , 0.15967134, 0.04300611,0.08314955, 0.03425328, 0.05417651, 0.04427107, 0.07115732, 0.04542704,0.16839289, 0.03772248, 0.08840031, 0.11684266, 0.0257566, 0.09902988])

}

#fittingを行う関数の読み出し