试图用pyhf组合一个教学示例

如何解决试图用pyhf组合一个教学示例

我试图学习有关pyhf的更多信息，而我对目标的理解可能会受到限制。我希望将我的HEP数据放在ROOT之外，但我可能对pyhf施加了期望，而这并不是作者打算使用的。

我想给自己写一个hello-world示例，但我可能只是不知道自己在做什么。我的误解也可能是我的统计知识方面的空白。

有了这个序言，让我解释一下我要探索的东西。

我有一些观察到的事件，可以为它们计算一些可观察的事件，并对该数据进行装箱的直方图。我假设存在两个贡献物理过程，我称之为信号和背景。我为这些过程生成了一些蒙特卡洛样本，从理论上讲，事件总数接近但不完全符合我的观察结果。

我想：

• 使数据适合这两个过程的假设
• 从适合中获得每个过程的事件数的最佳值
• 获取这些拟合值的不确定性
• 在适当的情况下，计算信号事件数的上限。

下面是我的入门代码，我所要做的只是适合ML，但我不确定要去哪里。我知道它并没有做我想做的事，但是我迷失了在RTD上找到的示例。我敢肯定是我，这不是对文档的批评。

import pyhf
import numpy as np
import matplotlib.pyplot as plt

nbins = 15

# Generate a background and signal MC sample`
MC_signal_events = np.random.normal(5,1.0,200)
MC_background_events = 10*np.random.random(1000)

signal_data = np.histogram(MC_signal_events,bins=nbins)[0]
bkg_data = np.histogram(MC_background_events,bins=nbins)[0]

# Generate an observed dataset with a slightly different
# number of events
signal_events = np.random.normal(5,180)
background_events = 10*np.random.random(1050)

observed_events = np.array(signal_events.tolist() + background_events.tolist())
observed_sample = np.histogram(observed_events,bins=nbins)[0]


# Plot these samples,if you like
plt.figure(figsize=(12,4))
plt.subplot(1,3,1)
plt.hist(observed_events,bins=nbins,label='Observations')
plt.legend()
plt.subplot(1,2)
plt.hist(MC_signal_events,label='MC signal')
plt.legend()
plt.subplot(1,3)
plt.hist(MC_background_events,label='MC background')
plt.legend()

# Use a very naive estimate of the background
# uncertainties
bkg_uncerts = np.sqrt(bkg_data)

print("Defining the PDF.......")
pdf = pyhf.simplemodels.hepdata_like(signal_data=signal_data.tolist(),\
                                     bkg_data=bkg_data.tolist(),\
                                     bkg_uncerts=bkg_uncerts.tolist())

print("Fit.......")
data = pyhf.tensorlib.astensor(observed_sample.tolist() + pdf.config.auxdata)

bestfit_pars,twice_nll = pyhf.infer.mle.fit(data,pdf,return_fitted_val=True)

print(bestfit_pars)
print(twice_nll)

plt.show()

解决方法

注意：此答案基于pyhf v0.5.2。

好的，因此看来您已经成功地确定了大部分重要零件。但是，有两种不同的方法可以执行此操作，具体取决于您喜欢的设置方式。在这两种情况下，我都假设您想要不受限制的适合度，并且您想要...

1. 使信号+背景模型适合观测数据

2. 使背景模型适合观察的数据

首先，让我们简要讨论不确定性。目前，张量背景默认为numpy，优化程序默认为scipy。查看文档：

• numpy backend
• scipy optimizer

但是，scipy优化器目前的一个不幸缺点是它无法返回不确定性。在适合之前（尽管我们通常建议尽早建议），您需要在代码中的任何地方进行操作，而不是使用minuit优化器：

pyhf.set_backend('numpy','minuit')

这将为您提供能够获得相关矩阵，拟合参数的不确定性以及粗麻布等漂亮功能。我们也正在努力使它与scipy保持一致，但是目前尚不成熟。

所有优化都通过我们的optimizer API进行，您目前可以通过文档中的mixin here查看。具体来说，签名是

minimize(
    objective,data,pdf,init_pars,par_bounds,fixed_vals=None,return_fitted_val=False,return_result_obj=False,do_grad=None,do_stitch=False,**kwargs)

这里有很多选择。让我们只关注一个事实，我们可以传递的关键字参数之一是return_uncertainties，它将通过为想要的拟合参数不确定性添加一列来更改最佳拟合参数。

---

1。信号+背景

在这种情况下，我们只想使用默认模型

result,twice_nll = pyhf.infer.mle.fit(
    data,return_uncertainties=True,return_fitted_val=True
)

bestfit_pars,errors = result.T

2。仅背景

在这种情况下，我们需要关闭信号。我们这样做的方法是将感兴趣的参数（POI）设置为0.0。然后，我们可以类似的方式获得仅背景模型的拟合参数，但使用fixed_poi_fit而不是无约束拟合：

result,twice_nll = pyhf.infer.mle.fixed_poi_fit(
    0.0,errors = result.T

请注意，这只是完成以下不受限制的拟合的一种简单方法

bkg_params = pdf.config.suggested_init()
fixed_params = pdf.config.suggested_fixed()

bkg_params[pdf.config.poi_index] = 0.0
fixed_params[pdf.config.poi_index] = True

result,init_pars=bkg_params,fixed_params=fixed_params,errors = result.T

---

希望这可以进一步澄清问题！

,

Giordon的解决方案应该回答您所有的问题，但是我想我也应该写出代码来基本上解决我们所能解决的所有问题。 我还可以自由更改一些值，以使信号不会太强，以至于观察到的CL值离巴西波段的右边并不远（结果并没有错，但是在那时谈论使用发现测试统计信息然后设置限制可能更有意义。：））

环境

在此示例中，我将设置一个干净的Python 3虚拟环境，然后安装依赖项（这里我们将使用pyhf v0.5.2）

$ python3 -m venv "${HOME}/.venvs/question"
$ . "${HOME}/.venvs/question/bin/activate"
(question) $ cat requirements.txt
pyhf[minuit,contrib]~=0.5.2
black
(question) $ python -m pip install -r requirements.txt

代码

虽然我们无法轻易获得信号事件数量和背景事件的最佳拟合值，但我们绝对可以推断出信号强度的最佳拟合值。 以下代码段（仅由于可视化而使之很长）应该解决您问题的所有要点。

# answer.py
import numpy as np
import pyhf
import matplotlib.pyplot as plt
import pyhf.contrib.viz.brazil

# Goals:
#    - Fit the model to the observed data
#    - Infer the best fit signal strength given the model
#    - Get the uncertainties on the best fit signal strength
#    - Calculate an 95% CL upper limit on the signal strength


def plot_hist(ax,bins,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.png")

    # 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:]}"
    )

    # 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.png")


if __name__ == "__main__":
    main()

运行时会给出

(question) $ python answer.py
best fit parameters:        
 * signal strength: 1.5884737977889158 +/- 0.7803435235862329        
 * nuisance parameters: [0.99020988 1.06040191 0.90488207 1.03531383 1.09093327 1.00942088
 1.07789316 1.01125627 1.06202964 0.95780043 0.94990993 1.04893286
 1.0560711  0.9758487  0.93692481 1.04683181 1.05785515 0.92381263
 0.93812855 0.96751869]        
 * nuisance parameter uncertainties: [0.06966439 0.07632218 0.0611428  0.07230328 0.07872258 0.06899675
 0.07472849 0.07403246 0.07613661 0.08606657 0.08002775 0.08655314
 0.07564512 0.07308117 0.06743479 0.07383134 0.07460864 0.06632003
 0.06683251 0.06270965]

Performing hypothesis tests

/home/stackoverflow/.venvs/question/lib/python3.7/site-packages/pyhf/infer/calculators.py:229: RuntimeWarning: invalid value encountered in double_scalars
  teststat = (qmu - qmu_A) / (2 * self.sqrtqmuA_v)
Observed Limit on µ: 2.89
-----
Expected (-2 σ) Limit on µ: 0.829
Expected (-1 σ) Limit on µ: 1.110
Expected        Limit on µ: 1.542
Expected (1 σ) Limit on µ: 2.147
Expected (2 σ) Limit on µ: 2.882

让我们知道您是否还有其他疑问！