Спасибо за хороший вопрос @ robsol90.
Если мы визуально проверяем содержимое модели (откройте файл ROOT и посмотрите на историограммы в TBrowser) или просто распечатайте содержимое ( после преобразования XML + ROOT в JSON)
>>> import json
>>> with open("meas.json") as spec_file:
... spec = json.load(spec_file)
...
>>> print(json.dumps(spec, indent=2, sort_keys=True))
мы видим, что в модели много корзин с нулями. Это проблема, так как HistFactory основана на пуассоновском алгоритме, и поскольку пуассоновское значение pmf определено строго для параметров скорости, превышающих 0, эти истинные 0 корзин вызовут ошибки (и они это делают). Однако, если мы просто проанализируем spe c и добавим очень небольшое смещение (epsilon), тогда подгонка может продолжаться без проблем. Таким образом, эта проблема на самом деле оказывается очень похожей на этот вопрос ( Ошибка сходимости Fit в pyhf для модели с малым сигналом ), и она не сразу становится очевидной.
Мы понимаем, что выбранная вами модель игрушки была Предполагается, что он минимален и прост, но поскольку в действительности вы почти никогда не встретите такой редкий регион анализа, эта проблема с игрушкой становится сложной. В будущем мы предпримем усилия, чтобы автоматически маскировать ячейки, являющиеся истинными нулями в модели, чтобы вообще избежать этой проблемы для пользователей.
Я также дам ниже некоторый код, который исправляет описанная выше проблема, а также некоторый дополнительный пример кода.
Во-первых, для ясности, давайте установим sh нашу среду
Среду
$ "$(which python3)" --version
Python 3.7.5
$ python3 -m venv "${HOME}/.venvs/question"
$ . "${HOME}/.venvs/question/bin/activate"
(question) $ cat requirements.txt
pyhf[xmlio]~=0.4.0
black
(question) $ python -m pip install -r requirements.txt
(question) $ root-config --version
6.18/04
Код
Давайте также разбить вещи на несколько этапов кода. Сначала давайте рассмотрим фрагмент кода от XML до ROOT, который я изменил, чтобы была более разумная выборка модели, отображаемой в наблюдаемых данных (, но не нужно было в качестве вашего оригинальный код будет работать и здесь).
# XML_to_ROOT.py
from ROOT import TH1F, TFile, TF1
def main():
left_bound = 95
right_bound = 115
n_bins = 80
# Model makeup
frac_bkg = 0.95
frac_sig = round(1.0 - frac_bkg, 2)
bkg_model = TF1("bkg_model", "TMath::Gaus(x,100,0.5,true)", left_bound, right_bound)
sig_model = TF1("sig_model", "TMath::Gaus(x,105,0.2,true)", left_bound, right_bound)
obs_model = TF1(
"obs_model",
f"({frac_bkg}*bkg_model)+({frac_sig}*sig_model)",
left_bound,
right_bound,
)
# Samples from model
n_sample = 10000
n_bkg = int(frac_bkg * n_sample)
n_sig = int(frac_sig * n_sample)
bkg_nominal = TH1F("bkg_nominal", "", n_bins, left_bound, right_bound)
bkg_nominal.FillRandom("bkg_model", n_bkg)
sig_nominal = TH1F("sig_nominal", "", n_bins, left_bound, right_bound)
sig_nominal.FillRandom("sig_model", n_sig)
data_nominal = TH1F("data_nominal", "", n_bins, left_bound, right_bound)
data_nominal.FillRandom("obs_model", n_sample)
meas = TFile("meas.root", "RECREATE")
bkg_nominal.Write()
sig_nominal.Write()
data_nominal.Write()
meas.Close()
if __name__ == "__main__":
main()
Теперь, чтобы упростить ситуацию позже, сгенерируем наш файл XML и ROOT, а затем преобразуем их в JSON spe c
(question) $ python XML_to_ROOT.py
(question) $ pyhf xml2json --output-file meas.json meas.xml
Теперь, наконец, давайте адаптируем код в вашем вопросе, чтобы убедиться, что ни одна ячейка в модели не содержит истинных 0 с, добавив все ячейки со смещением 1e-20 (просто чтобы продемонстрировать, что важно то, что они не равны нулю)
# answer.py
import os
import json
import pyhf.readxml
import numpy as np
def main():
with open("meas.json") as spec_file:
spec = json.load(spec_file)
# Pad true zeros to avoid error with evaluating Poisson(x|0)
epsilon = 1e-20
bkg = np.asarray(spec["channels"][0]["samples"][0]["data"]) + epsilon
sig = np.asarray(spec["channels"][0]["samples"][1]["data"]) + epsilon
spec["channels"][0]["samples"][0]["data"] = bkg.tolist()
spec["channels"][0]["samples"][1]["data"] = sig.tolist()
workspace = pyhf.Workspace(spec)
model = workspace.model(measurement_name="meas")
data = workspace.data(model)
best_fit_pars = pyhf.infer.mle.fit(data, model)
print(f"initialization parameters: {model.config.suggested_init()}")
print(
f"best fit parameters:\
\n * signal strength: {best_fit_pars[0]}\
\n * nuisance parameters: {best_fit_pars[1:]}"
)
if __name__ == "__main__":
main()
Теперь мы получаем
(question) $ python answer.py
initialization parameters: [1.0, 1.0, 1.0]
best fit parameters:
* signal strength: 1.000000316044688
* nuisance parameters: [0.99884051 1.02202245]
В качестве дополнительной демонстрации, что это действительно просто из-за истинных нулей, рассмотрите следующий пример с 2 бинами, который спроектирован так, чтобы завершиться с ошибкой.
# fail.py
import os
import json
import pyhf.readxml
import numpy as np
def main():
with open("meas.json") as spec_file:
spec = json.load(spec_file)
# Fails
bkg = np.asarray([0, 0])
sig = np.asarray([0, 1])
obs = np.asarray([1, 1])
# # Fails
# bkg = np.asarray([1, 0])
# sig = np.asarray([0, 0])
# obs = np.asarray([1, 1])
# # Fails
# bkg = np.asarray([0, 0])
# sig = np.asarray([0, 0])
# obs = np.asarray([1, 1])
# # Pass
# bkg = np.asarray([1e-9, 0])
# sig = np.asarray([0, 1e-9])
# obs = np.asarray([1, 1])
spec["channels"][0]["samples"][0]["data"] = bkg.tolist()
spec["channels"][0]["samples"][1]["data"] = sig.tolist()
spec["observations"][0]["data"] = obs.tolist()
workspace = pyhf.Workspace(spec)
model = workspace.model(measurement_name="meas")
data = workspace.data(model)
best_fit_pars = pyhf.infer.mle.fit(data, model)
if __name__ == "__main__":
main()