まずはCERN ROOT + C++ で実装する。
お手本通り、平均値0、標準偏差1、ガウス分布(正規分布)に沿う乱数を10000個作り、ROOTのヒストグラムに詰め、 TF1 の ガウス分布 gaus でフィッティングした。オプション等は何もつけていない。も普通に結果を引用した。
正規分布の関数は、
である。
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
ガウス分布の式は自分で作成し与えた。また、は
scipy.stats.chisquare を使って計算した。ここでポイントは、誤差を適切に与えることである。ビンあたりのエントリー数に対して、誤差
の配列を
curve_fit の sigma に与えるとよい。
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/
