フィッティングと桁落ち

スライドを作る上で、私の研究成果は今後どのような展開となるか、というのを示すために、実験データの統計が増えたと仮定した時の想定結果を出したくなり、久々に計算機にログイン。

私の研究では、仮想の統計での多次元フィッティング(約100のfloating parameter)で最終結果を出すが、何故かうまくいかないこと1~2時間。。。

素粒子実験の世界ではROOTという、CERNから出ている科学計算パッケージを使う習慣となっており、通常多次元フィッティングする場合はMinuitという最小値探索パッケージのお世話になる。Minuitは統計解析で流行りのマルコフ連鎖法は使わず、普通の勾配法を使う。これは単峰性が保証されていないことの方が多いからだと思われる。

それはさておき、色々ソースコードやログを眺めていたところ、correlation matrixの計算が不能になることでフィッティングが失敗に終わっているようだ。よく見ると、correlation matrixがnon positive definiteになるからエラーで終了しましたよ、と書いてある。Google検索してみると、あんまりこの手の情報でわかりやすくまとまっているものは無く、とりあえず原因はいろいろなケースがある、ということだけ分かった。あんまりこういうエラーに引っかかる人って少ないのだろうか。

いろいろ考えたが、私の場合統計を増やしてから失敗し始めたので、統計エラーが小さくなった分、より小さい桁での計算になった --> 恐らく桁落ちのせい、だと思って、ROOTのパッケージで計算機epsilonを大きめに設定したら幾つかのケースで上手く行き始めた。

ちなみに、計算機epsilonの設定は、

ROOT.Math.MinimizerOptions.SetDefaultTolerance( Scale )

という感じで行う。私の場合、fittingの後にp-valueが欲しいのでhypothesis testまでやるため、RooStats.AsymptoticCalculator を使うが、このソースコードとログを見比べる限り、セットする値は1.000000001とかの慣例的なものではなく、デフォルトの1.001に対するスケールを入れる。とりあえず私は10倍 (つまり1.01) にしてみたら、いい感じに計算不能だったfittingが成功し始めた。大体フィット結果も主要なパラメターは10倍する前とほぼ同じなので良しとする。ちなみにここでのtoleranceは - log likelihoodの最小値を出す時に、勾配の絶対値の値がその値を切るかどうか、ということだけに使われる。

統計を更に大幅に上げると、10倍のスケールでも上手く行かないので、少しずつ上げて試す必要がありそうだ。しかし１回のfittingでreal timeで30分くらいかかるので(CPU 20個並列に使った状態で)、何回も試すのは骨が折れそうだな、、、と思った。

それではこの辺で！