如何解决使用“ pyhf”浮动信号强度和背景强度
您好pyhf用户和开发人员！
我有一个来自previous question的问题，因此我将从响应之一和其次修改中提供的anwser.py代码开始。
因此，我使用响应中给定的参数运行拟合，但是然后我想查看拟合的结果，因此我添加了一些代码以根据拟合结果对原始模板进行加权，并用覆盖图重新绘制。完整的脚本在下面。
import pyhf
import pyhf.contrib.viz.brazil
import numpy as np
import matplotlib.pylab as plt
# - Get the uncertainties on the best fit signal strength
# - Calculate an 95% CL upper limit on the signal strength
tag = "ORIGINAL"
def plot_hist(ax,bins,data,bottom=0,color=None,label=None):
bin_width = bins[1] - bins[0]
bin_leftedges = bins[:-1]
bin_centers = [edge + bin_width / 2.0 for edge in bin_leftedges]
ax.bar(
bin_centers,bin_width,bottom=bottom,alpha=0.5,color=color,label=label
)
def plot_data(ax,label="Data"):
bin_width = bins[1] - bins[0]
bin_leftedges = bins[:-1]
bin_centers = [edge + bin_width / 2.0 for edge in bin_leftedges]
ax.scatter(bin_centers,color="black",label=label)
def invert_interval(test_mus,hypo_tests,test_size=0.05):
# This will be taken care of in v0.5.3
cls_obs = np.array([test[0] for test in hypo_tests]).flatten()
cls_exp = [
np.array([test[1][idx] for test in hypo_tests]).flatten() for idx in range(5)
]
crossing_test_stats = {"exp": [],"obs": None}
for cls_exp_sigma in cls_exp:
crossing_test_stats["exp"].append(
np.interp(
test_size,list(reversed(cls_exp_sigma)),list(reversed(test_mus))
)
)
crossing_test_stats["obs"] = np.interp(
test_size,list(reversed(cls_obs)),list(reversed(test_mus))
)
return crossing_test_stats
def main():
np.random.seed(0)
pyhf.set_backend("numpy","minuit")
observable_range = [0.0,10.0]
bin_width = 0.5
_bins = np.arange(observable_range[0],observable_range[1] + bin_width,bin_width)
n_bkg = 2000
n_signal = int(np.sqrt(n_bkg))
# Generate simulation
bkg_simulation = 10 * np.random.random(n_bkg)
signal_simulation = np.random.normal(5,1.0,n_signal)
bkg_sample,_ = np.histogram(bkg_simulation,bins=_bins)
signal_sample,_ = np.histogram(signal_simulation,bins=_bins)
# Generate observations
signal_events = np.random.normal(5,int(n_signal * 0.8))
bkg_events = 10 * np.random.random(int(n_bkg + np.sqrt(n_bkg)))
observed_events = np.array(signal_events.tolist() + bkg_events.tolist())
observed_sample,_ = np.histogram(observed_events,bins=_bins)
# Visualize the simulation and observations
fig,ax = plt.subplots()
fig.set_size_inches(7,5)
plot_hist(ax,_bins,bkg_sample,label="Background")
plot_hist(ax,signal_sample,bottom=bkg_sample,label="Signal")
plot_data(ax,observed_sample)
ax.legend(loc="best")
ax.set_ylim(top=np.max(observed_sample) * 1.4)
ax.set_xlabel("Observable")
ax.set_ylabel("Count")
fig.savefig("components_{0}.png".format(tag))
# Build the model
bkg_uncerts = np.sqrt(bkg_sample)
model = pyhf.simplemodels.hepdata_like(
signal_data=signal_sample.tolist(),bkg_data=bkg_sample.tolist(),bkg_uncerts=bkg_uncerts.tolist(),)
data = pyhf.tensorlib.astensor(observed_sample.tolist() + model.config.auxdata)
# Perform inference
fit_result = pyhf.infer.mle.fit(data,model,return_uncertainties=True)
bestfit_pars,par_uncerts = fit_result.T
print(
f"best fit parameters:\
\n * signal strength: {bestfit_pars[0]} +/- {par_uncerts[0]}\
\n * nuisance parameters: {bestfit_pars[1:]}\
\n * nuisance parameter uncertainties: {par_uncerts[1:]}"
)
# Visualize the results
fit_bkg_sample = []
for w,b in zip(bestfit_pars[1:],bkg_sample):
fit_bkg_sample.append(w*b)
fit_signal_sample = bestfit_pars[0]*np.array(signal_sample)
fig,fit_bkg_sample,fit_signal_sample,bottom=fit_bkg_sample,observed_sample)
ax.legend(loc="best")
ax.set_ylim(top=np.max(observed_sample) * 1.4)
ax.set_xlabel("Observable")
ax.set_ylabel("Count")
fig.savefig("components_after_fit_{0}.png".format(tag))
# Perform hypothesis test scan
_start = 0.0
_stop = 5
_step = 0.1
poi_tests = np.arange(_start,_stop + _step,_step)
print("\nPerforming hypothesis tests\n")
hypo_tests = [
pyhf.infer.hypotest(
mu_test,return_expected_set=True,return_test_statistics=True,qtilde=True,)
for mu_test in poi_tests
]
# Upper limits on signal strength
results = invert_interval(poi_tests,hypo_tests)
print(f"Observed Limit on µ: {results['obs']:.2f}")
print("-----")
for idx,n_sigma in enumerate(np.arange(-2,3)):
print(
"Expected {}Limit on µ: {:.3f}".format(
" " if n_sigma == 0 else "({} σ) ".format(n_sigma),results["exp"][idx],)
)
# Visualize the "Brazil band"
fig,5)
ax.set_title("Hypothesis Tests")
ax.set_ylabel(r"$\mathrm{CL}_{s}$")
ax.set_xlabel(r"$\mu$")
pyhf.contrib.viz.brazil.plot_results(ax,poi_tests,hypo_tests)
fig.savefig("brazil_band_{0}.png".format(tag))
if __name__ == "__main__":
main()
当我运行它时，得到以下图。第一个是原始的观测/模拟，第二个是根据拟合结果对模拟进行缩放。
所以这一切看起来都不错，我想我知道他发生了什么事。
但是现在我问一个问题，在您的示例中稍作改动会更好地说明这一问题。
我将修改生成的模拟和观察值，以使模拟和样本中存在不同数量的背景事件。我也在使信号更重要。这是一个示例，其中在进行拟合之前，我无法很好地估算出背景贡献。在您提供的示例中，模拟样本和数据的背景事件数相同，而在现实世界中情况并非如此。
因此，我转到上面的代码，并更改了这些行。
n_bkg = 2000
n_signal = 200
# Generate simulation
bkg_simulation = 10 * np.random.random(n_bkg)
signal_simulation = np.random.normal(5,int(n_signal * 0.8))
bkg_events = 10 * np.random.random(n_bkg - 300)
拟合度不是很好，并且我不希望这样，因为我锁定了背景事件的数量，对每个仓中的泊松波动进行了模运算。相关图（拟合之前/之后）显示在这里。
我可能曾想过另一种解决方法，那就是添加另一个无干扰的浮动参数，该参数表示背景强度，同时仍使各个垃圾箱在泊松波动范围内变化。出于这个原因，信号仓也不会（也不应该）波动吗？
在那种情况下，我将从模拟样本中大量的数据点开始，以获得更“真实”的分布（我知道这并不严格）。一旦拟合降低了信号/背景事件的数量，泊松波动将变得更加明显。
我敢肯定，似然函数的优化/最小化变得困难得多，但是如果我们锁定总体背景规范化，那感觉就像我们太早限制拟合。还是我想念什么？
一如既往地感谢您的帮助和答复！
解决方法
您应该可以通过向背景组件添加“ normfactor”修饰符来添加新的麻烦
例如
spec = {'channels': [{'name': 'singlechannel','samples': [{'name': 'signal','data': [10.0],'modifiers': [{'name': 'mu','type': 'normfactor','data': None}]},{'name': 'background','data': [50.0],'modifiers': [{'name': 'bkgnorm','data': None}]}]}]}
查看信号和背景之间的对称性
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
