我正在尝试构建一个非常小的pyhf示例：两个高斯，一个信号和一个背景，但我无法使其正常工作。我的python代码是：
import pyhf.readxml
import os
from ROOT import TH1F, TFile, TF1
mygaus = TF1("mygaus","TMath::Gaus(x,100,.5)",95, 115)
mygaus2 = TF1("mygaus2","TMath::Gaus(x,110,.2)",95, 115)
mygaus_data = TF1("mygaus_data","TMath::Gaus(x,110,.2)+TMath::Gaus(x,100,.5)",95, 115)
bkg_nominal = TH1F('bkg_nominal', '', 80, 95, 115)
bkg_nominal.FillRandom("mygaus", 10000)
sig_nominal = TH1F('sig_nominal', '', 80, 95, 115)
sig_nominal.FillRandom("mygaus2", 5000)
data_nominal = TH1F('data_nominal', '', 80, 95, 115)
data_nominal.FillRandom("mygaus_data", 10000)
meas = TFile('meas.root', 'RECREATE')
bkg_nominal.Write()
sig_nominal.Write()
data_nominal.Write()
meas.Close()
spec = pyhf.readxml.parse('meas.xml', os.getcwd())
workspace = pyhf.Workspace(spec)
pdf = workspace.model(measurement_name='meas')
data = workspace.data(pdf)
workspace.get_measurement(measurement_name='meas')
best_fit = pyhf.infer.mle.fit(data, pdf)
我基本上是从文档中的示例复制来的XML文件是这样写的
meas.xml
<!DOCTYPE Combination SYSTEM 'HistFactorySchema.dtd'>
<Combination OutputFilePrefix="workspace" >
<Input>./meas_channel1.xml</Input>
<Measurement Name="meas" Lumi='1' LumiRelErr='0.1' ExportOnly="False" >
<POI>signorm</POI>
</Measurement>
</Combination>
meas_channel1.xml
<!DOCTYPE Channel SYSTEM 'HistFactorySchema.dtd'>
<Channel Name="channel1" InputFile="" >
<Data HistoName="data_nominal" InputFile="meas.root" />
<StatErrorConfig RelErrorThreshold="0.05" ConstraintType="Gaussian" />
<Sample Name="bkg" HistoName="bkg_nominal" InputFile="meas.root" NormalizeByTheory="True" >
<NormFactor Name="bkgnorm" Val="1" High="3" Low="0" Const="False" />
</Sample>
<Sample Name="sig" HistoName="sig_nominal" InputFile="meas.root" NormalizeByTheory="True" >
<NormFactor Name="signorm" Val="1" High="3" Low="0" Const="False" />
</Sample>
</Channel>
它看起来非常简单，我能够绘制直方图。但是，当我收到此错误消息时：
ERROR:pyhf.optimize.opt_scipy: fun: nan
jac: array([nan, nan, nan])
message: 'Inequality constraints incompatible'
nfev: 5
nit: 1
njev: 1
status: 4
success: False
x: array([1., 1., 1.])
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
<ipython-input-14-54e7c2f0a645> in <module>
2 data = workspace.data(pdf)
3 workspace.get_measurement(measurement_name='meas')
----> 4 best_fit = pyhf.infer.mle.fit(data, pdf)
/usr/local/lib/python3.7/site-packages/pyhf/infer/mle.py in fit(data, pdf, init_pars, par_bounds, **kwargs)
34 init_pars = init_pars or pdf.config.suggested_init()
35 par_bounds = par_bounds or pdf.config.suggested_bounds()
---> 36 return opt.minimize(twice_nll, data, pdf, init_pars, par_bounds, **kwargs)
37
38
/usr/local/lib/python3.7/site-packages/pyhf/optimize/opt_scipy.py in minimize(self, objective, data, pdf, init_pars, par_bounds, fixed_vals, return_fitted_val)
45 )
46 try:
---> 47 assert result.success
48 except AssertionError:
49 log.error(result)
AssertionError:
这很奇怪，因为我没有任何不平等约束。我想我在做些蠢事，请您帮忙吗？谢谢！
参考方案
感谢您对@ 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基于Poisson以及Poisson p.m.f，这是一个问题。严格为速率参数定义大于0的值，这些真实的0 bin将导致错误（并且确实如此）。
但是，如果我们仅解析规格并添加一个非常小的偏移量（epsilon），则拟合就可以顺利进行。
因此，实际上这个问题最终变得非常类似于该问题（Fit convergence failure in pyhf for small signal model），而没有立即显现出来。
我们知道您设置的玩具模型应该是最小且容易的，但是实际上，您几乎永远不会遇到如此稀疏的分析区域，因此玩具问题变得很困难。尽管我们会在将来做出努力，以自动掩盖模型中为零的垃圾箱，以完全避免用户遇到此问题。
我还将在下面提供一些代码来解决您在上面遇到的问题，以及一些其他示例代码。
首先，非常清楚，让我们建立环境
环境
$ "$(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规范
(question) $ python XML_to_ROOT.py
(question) $ pyhf xml2json --output-file meas.json meas.xml
现在，最后，让我们修改问题代码，以确保所有模型中的bin都不包含真实的0，方法是将所有bin填充为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 bin示例，该示例经过设计可导致您的错误而失败。
# 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()
我正在使用python pandas处理一些数据。我已使用以下代码将数据导出到excel文件。writer = pd.ExcelWriter('Data.xlsx'); wrong_data.to_excel(writer,"Names which are wrong", index = False); writer.…
用大写字母拆分字符串，但忽略AAA Python Regex - python我的正则表达式：vendor = "MyNameIsJoe. I'mWorkerInAAAinc." ven = re.split(r'(?<=[a-z])[A-Z]|[A-Z](?=[a-z])', vendor) 以大写字母分割字符串，例如：'我的名字是乔。 I'mWorkerInAAAinc”变成…
Python-crontab模块 - python我正在尝试在Linux OS（CentOS 7）上使用Python-crontab模块我的配置文件如下：{ "ossConfigurationData": { "work1": [ { "cronInterval": "0 0 0 1 1 ?", "attribute&…
字符串文字中的正斜杠表现异常 - python为什么S1和S2在撇号位置方面表现不同？S1="1/282/03/10" S2="4/107/03/10" R1="".join({"N\'" ,S1,"\'" }) R2="".join({"N\'…
在Flask中测试文件上传 - python我在Flask集成测试中使用Flask-Testing。我有一个表单，该表单具有我要为其编写测试的徽标的文件上传，但是我不断收到错误消息：TypeError: 'str' does not support the buffer interface。我正在使用Python3。我找到的最接近的答案是this，但是它对我不起作用。这是我的许多尝…