python で fitting する簡単な例

概要

python でモデルフィットしてエラーを出すための簡単な方法を示す。極力短いコードで確実に動くサンプルを意図したので、エラー判定や細かな処理はしておりません。コードを見ればOKな人は私のgoogle Colab のページを参照ください。

フィットとは？

その前に、フィットとはどういう手続かを確認しておく。

データとモデルと用意する。ここでのモデルは、変数と定数を与えると値を返す関数である。
モデルの定数を与え、変数の初期値をセットする。
データとモデルの合具合を評価する関数を定義する。原理的には何でもよいが、カイ２乗値を使うことが多いので、この例ではそれに準じる。要するに、モデルとデータのズレの量である。初期値から計算を始め、ズレが最小になるように変数を初期値から動かし、ズレが極小になる場所を見つける(e.g, Levenberg-Marquardt Method)。ここではscipyの実装されたLM法 scipy.optimize.leastsq を用いる。(より一般的な場合は最後に補足)
個々のパラメータに付随する誤差および"モデルとの合い具合"を評価する。LM法では、エラーは変数の誤差の共分散行列の対角成分を信頼区間として評価するのでこれに準じた。注意としては、この方法では誤差は必ず正負で同じになり、科学論文で使えるエラーかどうかは考える必要ある。(e.g., X線衛星業界の話だと，xspec では error コマンドを使ってエラーを求めなさい、と同じ。)

これに準じてソースコードを書いた。

コードをみたらOKな人はこれを読んでください。

サンプルコード

直線フィットとガウシアンフィットの例を示す。

どちらもデータをベタにコードに書いて、初期値は適当にセットした。モデルの初期値は重要で、これが全然違うと収束しないことには注意する。この例ではそこそこな値を手打ちしてある。
モデルは mymodel という関数で設定する。下記の例では定数 consts は使ってないが、必要ならこれで定数を設定できる。パラメータの数や定数の数も任意に変更できるが、計算途中で発散が zero div が発生しないようには注意が必要。
mymodel関数の引数は xvalues,yvalues,params,constsの４つ用意しているが、yvaluesはこの例では使っていない。yの値にモデルが依存する場合に使う。
so.leastsq はカイ二乗値を計算するcalcchi関数、モデルの初期値param_init,calcchi関数の初期値以外のパラメータargs=(consts,model_func, xvalues, yvalues, yerrors)で必要である。またこの例は,yerrorsにゼロが含まれるときにzero div が発生してコケる。エラーがゼロにならないように使うこと。カウント数がゼロになるような場合はそもそもカイ二乗の最小化をすることが正しくない(詳細は後述を参考)。
エラーは誤差行列から出したエラーなので、科学的に〇〇%の信頼区間などを厳密に議論したい場合はこのエラーでよいかは検討が必要。中心値はLM法で求めたものがよく合っているなら、他の手法でも大して変わらないことが多い。

python での直線フィットの例

fit_linear.py

fit_linear.py


#!/usr/bin/env python

__version__= '1.0'
# This is created from the code made by M.Sawada during SXS tests. 

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'serif'
import numpy as np
import scipy as sp
import scipy.optimize as so

def calcchi(params,consts,model_func,xvalues,yvalues,yerrors):
    model = model_func(xvalues,yvalues,params,consts)
    chi = (yvalues - model) / yerrors
    return(chi)

# optimizer
def solve_leastsq(xvalues,yvalues,yerrors,param_init,consts,model_func):

    param_output = so.leastsq(
        calcchi,
        param_init,
        args=(consts,model_func, xvalues, yvalues, yerrors),
        full_output=True)
    param_result, covar_output, info, mesg, ier = param_output
    error_result = np.sqrt(covar_output.diagonal())
    dof = len(xvalues) - 1 - len(param_init)
    chi2 = np.sum(np.power(calcchi(param_result,consts,model_func,xvalues,yvalues,yerrors),2.0))
    return([param_result, error_result, chi2, dof])

# model
def mymodel(xvalues,yvalues,params,consts): 
    p0,p1 = params
    c0,c1 = consts
    model_y = p0*xvalues + p1
#    model_y = p0*np.exp(-0.5*((xvalues-p1)/p2)**2.) + c1
#    model_y = p0*np.exp(-0.5*((xvalues-p1)/p2)**2.) + c0*xvalues + c1
    return(model_y)
            
# set data
x  = np.array([ 1, 2, 3,  4,  5, 6, 7,  8,  9, 10, 11, 12, 13, 14, 15])
y  = np.array([ 1, 2, 5,  4,  8, 5, 9, 10, 13, 11, 15,  9, 10, 17, 25])
ye = np.array([ 3, 3, 2,  2,  3, 2, 3,  5,  6, 1,  3,  1,  3,  3,  2])

# initialize parameter 
a_init = 1.
b_init = 0.5
init_params = np.array([a_init, b_init])
# initialize constants 
consts = np.array([0,0]) 

# do fit 
result, error, chi2, dof = solve_leastsq( x, y, ye, init_params, consts, mymodel)

# get results 
a = result[0]
b = result[1]
ae = error[0]
be = error[1]

label = "Y = " + str("%4.2f(+/-%4.2f) " % (a,ae) ) + " X + " + str("%4.2f(+/-%4.2f)" % (b,be) )
print(label)

# plot results
plt.title("Linear fit : " + label)
model_y = mymodel(x,y,result,consts)
plt.errorbar(x, y, yerr = ye, fmt="bo")
plt.plot(x, model_y, 'r-')

plt.savefig("linearfit.png")
plt.show()

結果の絵がこちら。

python を用いたガウシアンフィットの例

fit_gauss.py

fit_gauss.py


#!/usr/bin/env python

__version__= '1.0'
# This is created from the code made by M.Sawada during SXS tests. 

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'serif'
import numpy as np
import scipy as sp
import scipy.optimize as so

def calcchi(params,consts,model_func,xvalues,yvalues,yerrors):
    model = model_func(xvalues,yvalues,params,consts)
    chi = (yvalues - model) / yerrors
    return(chi)

# optimizer
def solve_leastsq(xvalues,yvalues,yerrors,param_init,consts,model_func):

    param_output = so.leastsq(
        calcchi,
        param_init,
        args=(consts,model_func, xvalues, yvalues, yerrors),
        full_output=True)
    param_result, covar_output, info, mesg, ier = param_output
    error_result = np.sqrt(covar_output.diagonal())
    dof = len(xvalues) - 1 - len(param_init)
    chi2 = np.sum(np.power(calcchi(param_result,consts,model_func,xvalues,yvalues,yerrors),2.0))
    return([param_result, error_result, chi2, dof])

# model
def mymodel(xvalues,yvalues,params,consts): 
    p0,p1,p2 = params
    c0,c1 = consts
#    model_y = p0*xvalues + p1 
    model_y = p0*np.exp(-0.5*((xvalues-p1)/p2)**2.) + c1
#    model_y = p0*np.exp(-0.5*((xvalues-p1)/p2)**2.) + c0*xvalues + c1
    return(model_y)

# set data
x = np.array([ 1, 2, 3,  4,  5, 6,   7,  8,  9, 10, 11, 12, 13, 14, 15])
y = np.array([ 1, 2, 5, 10, 28, 32, 37, 51, 41, 25,  8,  3,  1,  2,  1])

# initialize parameter 
sigma_init = 5.
norm_init = 10.
init_params = np.array([norm_init, np.mean(x), sigma_init ])
# initialize constants 
consts = np.array([0,0]) 

# do fit 
result, error, chi2, dof = solve_leastsq( x, y, np.sqrt(y), init_params, consts, mymodel)

# get results 
ene = result[1]
sigma = np.abs(result[2])

enee = error[1]
sigmae = np.abs(error[2])

fwhm = 2.35 * sigma
fwhme = 2.35 * sigmae

label = "E = " + str("%4.2f(+/-%4.2f)" % (ene,enee)) + " FWHM = " + str("%4.2f(+/-%4.2f)" % (fwhm,fwhme) + " (FWHM)")
print(label)

# plot results
plt.title("Gauss fit : " + label)
model_y = mymodel(x,y,result,consts)
plt.errorbar(x, y, yerr = np.sqrt(y), fmt="bo")
plt.plot(x, model_y, 'r-')

plt.savefig("gaussfit.png")
plt.show()

結果の絵がこちら。あるいは、google Colab のページを参照ください。

モデルがガタガタしてみえるのは、簡単のため、データ点と同じ数だけモデルも描画しているため。例えば、numpy.interp などでデータ点数を増やしてモデルを表示すると綺麗になる。使い方は，y = np.interp(サンプリングしたいX座標の配列, 既知のX座標の配列, 既知のY座標の配列) のように，既知のx, y のペアを入れて，知りたい x 座標を入れると，知りたい y の値がかえってくる．これは単なる台形補完による内層なので，予想された範囲内の補完値しかでてこないので実験等では重宝します．

エラーについての補足

誤差は、numpy.diagonal で誤差の分散共分散行列の対角成分のルートをとってつけている。これは、エラーがガウス分布すると仮定すると、カイ２乗分布するパラメータの Hessian 行列の逆行列が、パラメータが真値(極小値)に近い時に、パラメータの誤差の分散共分散行列に一致するためである。

http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/DAVIES1/rd_bhatt_cvonline/node9.html
https://www.psychologie.uni-heidelberg.de/ae/meth/team/mertens/blog/hessian.nb.html
Hessian and Covariance Matrix 「Hessian matrix is equal to the inverse of the covariance matrix」 Hessian 行列の対角成分は、誤差行列(ガウシアンのシグマの２乗が対角に詰まった)の逆数になるので、Hessian 行列の逆行列のルートをとると標準偏差(~=パラメータの誤差)が得られる。

Hessian の逆行列(= sigma^2) 生成 ==> covariant Matrix を生成 ==> covariant Matrix の対角成分が誤差行列となる。 ２つの異なるパラメータの誤差や相関などは一切考慮してないことに注意

別の言い方としては、Fisher 情報量と Cramér–Rao (クラメール-ラオ)の不等式とか考えると、

何か下限値があって、それは確率分布関数の極値近傍の曲率に相当するような量(の逆数）が適切なような気がするので、+/-で対照的な誤差をつけるということは、曲率でえいやっと誤差をつけたことに相当するのだと思ってます。

より一般的な場合についての補足

この例は最も簡単なフィットのケースである。例えば、ポアソン統計が必要な場合は、ガウス分布ではなくてポアソン分布を使う必要があるので、カイ２乗の最小化問題ではなくなる。これにより、より一般的な非線形関数の最小化問題になる(カイ２乗の場合は二乗する関数の最小化という事前情報を使えたため高速で計算できた)。興味のある方は、下記のいつくかの文献を参考にしてほしい。

xspec の C 統計 : https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html　
ガウス統計とポアソン統計の比較の例 : https://academic.oup.com/pasj/article-abstract/71/4/75/5536242
計算方法 : https://link.springer.com/article/10.1007/s10909-014-1098-4

そのほか、微分不可能な関数の最小化が必要な場合は、

など参考にしてください。

python で直線フィットとガウシアンフィットをする簡単な方法