物理の駅 Physics station by 現役研究者

テクノロジーは共有されてこそ栄える

トップ > Python > CERN ROOTとPython scipy、lmfit を使って正規分布でフィッティングする方法を比較

CERN ROOTとPython scipy、lmfit を使って正規分布でフィッティングする方法を比較

Python C/C++ Fitting

同じデータを使って、CERN ROOT、Pythonのscipy、lmfit を使って正規分布でフィッティングするコードと結果を比較する。

まずは、データ生成部分

C++用

    std::vector<double> vx{-3.65,-3.45,-3.35,-3.25,-3.15,-3.05,-2.95,-2.85,-2.75,-2.65,-2.55,-2.45,-2.35,-2.25,-2.15,-2.05,-1.95,-1.85,-1.75,-1.65,-1.55,-1.45,-1.35,-1.25,-1.15,-1.05,-0.95,-0.85,-0.75,-0.65,-0.55,-0.45,-0.35,-0.25,-0.15,-0.05,0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95,1.05,1.15,1.25,1.35,1.45,1.55,1.65,1.75,1.85,1.95,2.05,2.15,2.25,2.35,2.45,2.55,2.65,2.75,2.85,2.95,3.05,3.15,3.25,3.35,3.55,3.75};
    std::vector<double> vy{1,1,2,3,4,6,1,7,15,8,18,22,20,40,36,41,68,75,92,117,130,118,151,185,192,205,297,280,302,330,314,362,388,366,403,388,423,367,416,389,358,338,329,296,294,283,236,208,171,153,129,116,91,65,81,52,41,33,31,18,23,23,14,8,9,6,2,2,1,3,1,2};
    std::vector<double> vyErrors;
    for (int i=0; i<vx.size(); i++){
        vyErrors.push_back(std::sqrt(vy[i]));
    }

Python用

vx = [-3.65,-3.45,-3.35,-3.25,-3.15,-3.05,-2.95,-2.85,-2.75,-2.65,-2.55,-2.45,-2.35,-2.25,-2.15,-2.05,-1.95,-1.85,-1.75,-1.65,-1.55,-1.45,-1.35,-1.25,-1.15,-1.05,-0.95,-0.85,-0.75,-0.65,-0.55,-0.45,-0.35,-0.25,-0.15,-0.05,0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95,1.05,1.15,1.25,1.35,1.45,1.55,1.65,1.75,1.85,1.95,2.05,2.15,2.25,2.35,2.45,2.55,2.65,2.75,2.85,2.95,3.05,3.15,3.25,3.35,3.55,3.75]
vy = [1,1,2,3,4,6,1,7,15,8,18,22,20,40,36,41,68,75,92,117,130,118,151,185,192,205,297,280,302,330,314,362,388,366,403,388,423,367,416,389,358,338,329,296,294,283,236,208,171,153,129,116,91,65,81,52,41,33,31,18,23,23,14,8,9,6,2,2,1,3,1,2]
vyError=[y**0.5 for y in vy]

次に、フィッティングの部分

ROOT CERN のコード。通常はヒストグラムに対してFitすることが多いが、今回はTGraphErrorsクラスを使う。

void FitGaussian() {
    //C++用のデータ生成部分をコピペ
    TGraph* graph = new TGraphErrors(vx.size(), vx.data(), vy.data(), nullptr, vyErrors.data());
    TF1* func = new TF1("func", "gaus", -5, 5);
    func->SetParameter(0, 400); // Constantの初期値
    func->SetParameter(1, 0); // Meanの初期値
    func->SetParameter(2, 1); // Sigmaの初期値
    graph->Fit("func", "Q");

    cout << "Constant " << func->GetParameter(0);
    cout << " +/- " << func->GetParError(0) << endl;
    cout << "Mean     " << func->GetParameter(1);
    cout << " +/- " << func->GetParError(1) << endl;
    cout << "Sigma    " << func->GetParameter(2);
    cout << " +/- " << func->GetParError(2) << endl;
    cout << "Chi2     " << func->GetChisquare() << endl;
    cout << "ndf      " << func->GetNDF() << endl;    
}

Python scipy.optimize.curve_fit を使うコード。パラメータの誤差を正しく計算するためabsolute_sigma =True とする。

import numpy as np
import pandas as pd
from scipy.optimize import curve_fit
from scipy.stats import chisquare

#Python用のデータ生成部分をコピペ
    
def gaussian_func(x, constant, mean, sigma):
    return constant * np.exp(- (x - mean) ** 2 / (2 * sigma ** 2))
parameter_initial = np.array([400, 0, 1]) #初期値を与える
popt, pcov = curve_fit(gaussian_func, vx, vy, sigma=vyError, absolute_sigma =True, p0=parameter_initial) #フィッティング
vy_fit = gaussian_func(np.array(vx), popt[0], popt[1], popt[2]) #最確値を使ってvyを計算
chisq, p_value = chisquare(vy_fit, f_exp=vy/np.mean(vy)*np.mean(vy_fit), ddof = 2) #カイ二乗を計算

stderr = np.sqrt(np.diag(pcov)) #分散共分散行列の対角成分からパラメータの誤差を計算
mat = np.vstack((popt,stderr)).T #最確値と誤差の行列
print("Constant {:.3f} +/- {:.5f}".format(*mat[0]))
print("Mean     {:.8f} +/- {:.8f}".format(*mat[1]))
print("Sigma    {:.6f} +/- {:.8f}".format(*mat[2]))
print("Chi2     {:.4f}".format(chisq))

Python lmfit ライブラリを使うコード。パラメータの誤差を正しく計算するため、重み weights には誤差の逆数を与えて scale_covar=Falseとする。

import numpy as np
import lmfit

#Python用のデータ生成部分をコピペ

model = lmfit.models.GaussianModel() #モデルの定義
params = model.guess(np.array(vy), x=np.array(vx)) #初期値を推定
result = model.fit(vy, params, x=vx, weights=1.0/np.array(vyError), scale_covar=False) #フィッティング

params = result.result.params
print("Constant {:.3f} +/- {:.5f}".format(params["height"].value,params["height"].stderr))
print("Mean     {:.8f} +/- {:.8f}".format(params["center"].value,params["center"].stderr))
print("Sigma    {:.6f} +/- {:.8f}".format(params["sigma"].value,params["sigma"].stderr))
print("Chi2     {:.4f}".format(result.chisqr))
print("ndf      {}".format(result.nfree))

各コードで計算した結果を比較すると、ほぼ一致している。

ROOT の CERNで fitting の結果

Constant 402.531 +/- 4.95816
Mean     -0.00397889 +/- 0.00989405
Sigma    0.983115 +/- 0.00704207
Chi2     83.6418
ndf      69

Python scipy.optimize.curve_fit の結果

Constant 402.531 +/- 4.96311
Mean     -0.00397894 +/- 0.00989406
Sigma    0.983115 +/- 0.00706389
Chi2     83.6418

Python lmfit ライブラリを使った結果

Constant 402.531 +/- 4.96310
Mean     -0.00397896 +/- 0.00989405
Sigma    0.983115 +/- 0.00706388
Chi2     83.6418
ndf      69

どのコードを使っても結果は変わらないので好きなコードを使うと良い。

lmfitはビルトイン(定義済み)関数が多く、カイ二乗なども勝手に計算してくれるので便利である。複雑な非線形関数を使ってfittingしたいときは scipyの curve_fit の方が便利だろう。