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

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

Python scipyのcurve_fit で導出したパラメータと標準誤差をCERN ROOTの結果と比較

まずはCERN ROOT + C++ で実装する。

お手本通り、平均値0、標準偏差1、ガウス分布(正規分布)に沿う乱数を10000個作り、ROOTのヒストグラムに詰め、 TF1 の ガウス分布 gaus でフィッティングした。オプション等は何もつけていない。\chi^{2}も普通に結果を引用した。

正規分布の関数は、


f(x)=const. \exp \! \left( -\frac{(x-mean)^2}{sigma^2} \right)

である。

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

次に、全く同じデータを使うために出力しておいたテキストデータを使って、 scipy.optimize.curve_fit でフィッティングさせる。 テキストデータはC++で発生させることが可能だが、C++のコンパイル環境がない人のためにファイルも提供する。

ガウス分布のフィッティング用のサンプルデータ Sample data

ガウス分布の式は自分で作成し与えた。また、\chi^{2}scipy.stats.chisquare を使って計算した。ここでポイントは、誤差を適切に与えることである。ビンあたりのエントリー数Nに対して、誤差\sqrt{N}の配列を curve_fitsigma に与えるとよい。

Chi2 83.641822
            Estimate  Std. error
Constant  402.530641    4.963109
Mean       -0.003979    0.009894
Sigma       0.983115    0.007064

推定値は完全に一致。誤差は若干異なるがなぜだろう。愛嬌でしょうか。誤差の3桁目までは一致しているので、実用上は問題ないだろう。 scipy.optimize.curve_fit においても、誤差を適切に与えることでROOTの結果と推定値は完全に一致することが分かった。

p値(p-value)を求めるときは、カイ二乗をc、自由度数をndfとすると、ROOTだとTMath::Prob(c, ndf)、Pythonだと

from scipy.stats import chi2
chi2.sf(c, ndf)

とする。今回の例だと TMath::Prob(83.641822, 69)=chi2.sf(83.641822, 69)=0.11055024 である。

以下ソースコード。

C++ と ROOT版 C++ と ROOT5版 ガウス分布をフィッティング ($1781181) · Snippets · GitLab

Python3版 Python3版 ガウス分布をフィッティング ($1781185) · Snippets · GitLab

ウェブの情報をいっぱい参考にしたが、その一部。

https://omedstu.jimdo.com/2018/07/16/python%E3%81%A7gaussian-fitting/

scipy.stats.chisquare — SciPy v1.11.1 Manual