物理の駅 by onsanai

Physics-station 研究で日々感じたことを忘れないための備忘録

Python3でROOT+C++と同様にフィッティングとパラメータの標準誤差を算出する

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

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

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

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

https://gitlab.com/snippets/1781539/raw?inline=false

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

f:id:onsanai:20181121122953p:plain:w500

Chi2 83.641822
            Estimate  Std. error
Constant  402.530641    5.464381
Mean       -0.003979    0.010893
Sigma       0.983115    0.007777

推定値は完全に一致。誤差は若干異なるがなぜだろう。愛嬌でしょうか。とにかく、 scipy.optimize.curve_fit においても、誤差を適切に与えることでROOTの結果と推定値は完全に一致することが分かった。以下ソースコード

C++ と ROOT版 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.1.0 Reference Guide