<jupyter_start><jupyter_text>My First Likelihood🎶 I'm the very Model of a simple HEP-like measurement... 🎶 HEP-like?So what do we do as experimentalists in High Energy Physics? We have a gigantic detector that serves as a counting machine for physics. We smash particles and take pictures of the aftermath. Then we propose hypotheses for new physics that we want to test using this data. So we will have, at the very least:* signal* background (with some uncertainty)* some observationsThis is probably the simplest model one could use in HistFactory... so let's make it.<jupyter_code>import json
import numpy as np
import matplotlib.pyplot as plt

import pyhf
from pyhf.contrib.viz import brazil  # not imported by default!
model = pyhf.simplemodels.hepdata_like(
    signal_data=[5.0, 10.0], bkg_data=[50.0, 60.0], bkg_uncerts=[5.0, 12.0]
)
model<jupyter_output><empty_output><jupyter_text>What did we just make? This returns a [`pyhf.Model`](https://pyhf.readthedocs.io/en/v0.6.0/_generated/pyhf.pdf.Model.htmlpyhf.pdf.Model) object. Let's check out the specification.<jupyter_code>print(json.dumps(model.spec, indent=2))<jupyter_output><empty_output><jupyter_text>As we can see from the specification, we defined a single two-bin channel ('singlechannel') with two samples: a **signal sample** and a **background sample**. The signal sample has an _unconstrained normalization factor_ $\mu$, the signal strength, while the background sample carries an _uncorrelated shape systematic_ controlled by parameters $\gamma_{1}$ and $\gamma_{2}$. The _background uncertainty_ for the bins is $10\%$ and $20\%$ respectively.These uncertainties are **absolute** (not relative!). We have 10% and 20% relative uncertainty for the bins in the background sample respectively.(*Note: we have a workspace defined with this simple model in `data/2-bin_1-channel.json` which we will look at later.*)Let's explore the model a little bit.<jupyter_code>print(f"  channels: {model.config.channels}")
print(f"     nbins: {model.config.channel_nbins}")
print(f"   samples: {model.config.samples}")
print(f" modifiers: {model.config.modifiers}")
print(f"parameters: {model.config.parameters}")
print(f"  nauxdata: {model.config.nauxdata}")
print(f"   auxdata: {model.config.auxdata}")<jupyter_output><empty_output><jupyter_text>Auxiliary Data Whoa, hold up! What's with the auxiliary data? Recall the HistFactory definition$$p(n,a|\theta) = \underbrace{\prod_{\mathrm{channel}\ c}\prod_{\mathrm{bin}\ b} \mathrm{Pois}(n_{cb} | \lambda_{cb}(\theta))}_{\text{main}} \underbrace{\prod_{\mathrm{constraint}\ \chi} p_\chi(a_\chi | \chi)}_{\text{auxiliary}} \qquad \lambda_{cb}(\theta) = \sum_{\mathrm{sample}\ s} \lambda_{cbs}(\theta)$$and the auxiliary data here is passed into a constraint function? There's two things going on here:1. the uncorrelated shapesys modifier is constrained by a Poisson2. the auxiliary data is fully determined by the shapesys 'data' and the 'background' dataSo if we were to explicitly write out the likelihood we were seeing symbolically, it would be:$$\begin{align}p_\text{main} \cdot p_\text{aux} &= \text{Pois}(n | \lambda) \cdot p_\text{aux}\\                                 &= \text{Pois}(n | \mu s + \gamma b) \cdot p_\text{aux}\\                                 &= \text{Pois}(n | \mu s + \gamma b) \cdot \text{Pois}(n_\text{aux} = (b/\delta b)^2 | \mu = (b/\delta b)^2 \gamma )\\\end{align}$$where $n = \{n_1, n_2\}$ for a 2-bin model (we're being slightly fast and loose with our mathematical notation here), and similarly for $s$, $b$, and $\gamma$.The 'shapesys' is defined in the [HistFactory paper](https://cds.cern.ch/record/1456844)... however, it can be a little hard to extract out the necessary information. We've provided a nice table of [Modifiers and Constraints](https://pyhf.readthedocs.io/en/v0.6.0/intro.htmlid24) in the introduction of our pyhf documentation to use as reference.![modifiers and constraints](assets/modifiers_and_constraints.png) Let's look at the first row here which is our uncorrelated shape modifier. This is a multiplicative modifier denoted by $\kappa$ per-bin (denoted by $\gamma_b$). Notice that the input for the constraint term requires $\sigma_b$ which is the relative uncertainty of that modifier. This is Poisson-constrained by $\sigma_b^{-2}$. Let's quickly verify by hand to convince ourselves of what's going on here:<jupyter_code>(np.array([5.0, 12.0]) / np.array([50.0, 60.0])) ** -2
model.config.auxdata<jupyter_output><empty_output><jupyter_text>Whew. Ok, I guess that makes sense. Data and ParametersWhat about the actual data for the model then?<jupyter_code>model.expected_data()<jupyter_output><empty_output><jupyter_text>Hrmm, that didn't work. It seems we need to specify the parameter values to get the data for the model at those parameter values. Well, we know the default for $\mu=1$ and $\gamma = 1$. Also recall that the uncorrelated shape allocates a parameter per bin, so our model has **3 parameters** with **2 modifiers**. So let's just try that:<jupyter_code>model.expected_data([1.0, 1.0, 1.0])<jupyter_output><empty_output><jupyter_text>This returns the data for the entire likelihood, the main model as well as the constraint (or auxiliary) model. We can also drop the auxdata to get the actual data.<jupyter_code>model.expected_data([1.0, 1.0, 1.0], include_auxdata=False)<jupyter_output><empty_output><jupyter_text>Similarly, there are also methods separately for the actual data and the auxdata, both which take in all the parameters:<jupyter_code>model.expected_actualdata([1.0, 1.0, 1.0])
model.expected_auxdata([1.0, 1.0, 1.0])<jupyter_output><empty_output><jupyter_text>The Parameter Ordering So how do we know what these parameters correspond to? The way `pyhf` works under the hood, you can imagine an N-dimensional tensor that is purely bookkeeping all the different expected rates and modifications (read: systematics) involved that we can apply tensor algebra calculations to. This means that the ordering of the entries in each dimension is important, so we need to keep track of the order of the channels, the samples, the parameters. In `pyhf`, you can ask for the parameters or modifiers in a model:<jupyter_code>model.config.parameters
model.config.modifiers<jupyter_output><empty_output><jupyter_text>However, these might not necessarily correspond with the order of the parameters when we build this tensor internally, and we keep this order under a different variable<jupyter_code>model.config.par_order<jupyter_output><empty_output><jupyter_text>The Parameter Multiplicity But this clearly doesn't answer everything, because we were using three numbers assigned to two parameters! How does that work? Well, recall that the only channel in this model has 2 bins:<jupyter_code>model.config.channel_nbins<jupyter_output><empty_output><jupyter_text>and one of the parameters (the uncorrelated shape systematic) will be configured by $\gamma_b$ - that is, one parameter for each bin of the sample in that channel. How does `pyhf` know that the parameter set is associated with two parameters? We can ask for information about the `paramset`:<jupyter_code>model.config.param_set("uncorr_bkguncrt")<jupyter_output><empty_output><jupyter_text>which we can see is constrained by a Poisson and<jupyter_code>model.config.param_set("uncorr_bkguncrt").n_parameters<jupyter_output><empty_output><jupyter_text>The Parameter Defaults So how do we know what we should set the parameters to? The nice thing is that this is up to you! The model will come with some suggested defaults for parameters...<jupyter_code>model.config.suggested_init()<jupyter_output><empty_output><jupyter_text>... and their bounds...<jupyter_code>model.config.suggested_bounds()<jupyter_output><empty_output><jupyter_text>... and if the parameter should be held constant in a fit or not ...<jupyter_code>model.config.suggested_fixed()<jupyter_output><empty_output><jupyter_text>So armed with this knowledge, we could have also done<jupyter_code>init_pars = model.config.suggested_init()
model.expected_actualdata(init_pars)<jupyter_output><empty_output><jupyter_text>To get the expected actual data from the model (signal + background). We could turn off the signal, corresponding with the normalization factor. We just need to know the index of the parameter:<jupyter_code>model.config.poi_index<jupyter_output><empty_output><jupyter_text>and then just change the parameter<jupyter_code>bkg_pars = [*init_pars]  # my clever way to "copy" the list by value
bkg_pars[model.config.poi_index] = 0
model.expected_actualdata(bkg_pars)<jupyter_output><empty_output><jupyter_text>Simple InferenceThe core of statistical analysis is the statistical model. For inference, it's viewed as a function of the parameters with the data fixed.$$\log L(\theta | x) = \log p(x | \theta)$$The value of the likelihood is a float. Let's try it for both the background-only model as well as the signal+background model.<jupyter_code>observations = [52.5, 65.0] + model.config.auxdata  # this is a common pattern!

model.logpdf(pars=bkg_pars, data=observations)
model.logpdf(pars=init_pars, data=observations)<jupyter_output><empty_output><jupyter_text>We're not performing inference just yet. We're simply computing the 'logpdf' of the model specified by the parameters $\theta$ against the provided data. To perform a fit, we use the [inference API](https://pyhf.readthedocs.io/en/v0.6.0/api.htmlinference) via `pyhf.infer`.To fit a model to data, we usually want to find the $\hat{\theta}$ which refers to the "Maximum Likelihood Estimate". This is often referred to mathematically by$$\hat{\theta}_\text{MLE} = \text{argmax}_\theta L(\theta | x)$$Let's perform a fit on this model to the provided observations we just made up.<jupyter_code>pyhf.infer.mle.fit(data=observations, pdf=model)<jupyter_output><empty_output><jupyter_text>So what can we say? With nominal signal `[5, 10]` and nominal background = `[50, 60]`, an observed count of `[52.5, 65]` suggests best fit values:* $\hat{\mu} \approx 0.5$,* $\hat{\gamma} \approx [1,1]$.  Simple Hypothesis TestingGreat, so we can fit. But that's not usually the interesting part, since we cannot make any statements about this. What we prefer to do is perform a hypothesis test to either:* reject Standard Model hypothesis - a discovery fit* reject Beyond the Standard Model hypothesis - an exclusion fitSo we need to compute test statistics in order to evaluate the hypothesis. If we know the distribution of the test statistic under two different hypotheses, we can compute a [CLs](https://en.wikipedia.org/wiki/CLs_method_(particle_physics)) value. This is a modified pseudo-frequentist *p*-value.$$\text{CL}_\text{s} = \frac{\text{CL}_\text{s+b}}{\text{CL}_\text{b}}$$with$$\text{CL}_\text{s+b} = \int_{t_\text{obs}}^\infty \text{d}t\ p_0(t | H_0) \qquad\qquad \text{CL}_\text{b} = \int_{t_\text{obs}}^\infty \text{d}t\ p_1(t | H_1)$$This is done by splitting the model parameters up into two groups:* parameters of interest (POI - can have multiple!)* nuisance parameters (NP)Recall that additionally, all parameters are either constrained or unconstrained. Parameters of interest happen to always be unconstrained parameters.$$    f(\boldsymbol{x}|\phi) = f(\boldsymbol{x}|\overbrace{\eta}^{\llap{\text{free}}},\underbrace{\chi}_{\llap{\text{constrained}}}) = f(\boldsymbol{x}|\overbrace{\psi}^{\rlap{\text{parameters of interest}}},\underbrace{\theta}_{\rlap{\text{nuisance parameters}}})$$In the vast majority of cases, the test statistic is the profile likelihood ratio, or variations of it:$$t_\psi(x) = \frac{L(\psi, \hat{\hat{\theta}})}{L({\hat{\psi}, \hat{\theta}})}$$where* $\hat{\hat{\theta}}$ is the best fitted value of the nuisance parameters (for fixed POIs)* $\hat{\psi}$ and $\hat{\theta}$ are the best fitted values in a global fitSo let's run a hypothesis test for* null hypothesis ($\mu = 1$) - "SUSY is real"* alternate hypothesis ($\mu = 0$) - "Standard Model explains it all"<jupyter_code>CLs_obs, CLs_exp = pyhf.infer.hypotest(
    1.0,  # null hypothesis
    [52.5, 65.0] + model.config.auxdata,
    model,
    test_stat="q",
    return_expected_set=True,
)
print(f"      Observed CLs: {CLs_obs:.4f}")
for expected_value, n_sigma in zip(CLs_exp, np.arange(-2, 3)):
    print(f"Expected CLs({n_sigma:2d} σ): {expected_value:.4f}")<jupyter_output><empty_output><jupyter_text>Nicely enough, `pyhf` is smart enough to let us know that, if for a model where $\hat{\psi} < 0$, we're not going to be able to catch it because our default bounds were bounding the parameter of interest from below at 0. For more information, see Equation 14 in [[1007.1727]](https://arxiv.org/abs/1007.1727).<jupyter_code>model.config.suggested_bounds()[model.config.poi_index]<jupyter_output><empty_output><jupyter_text>We should heed the warning and switch our test statistic in this case to a variation of it $q_\mu \to \tilde{q}_\mu$. We could also change the bounds from the default to allow $\mu$ to float below zero. In this case, we're just going to use a more appropriate test statistic as we know $\hat{\mu} \geq 0$.<jupyter_code>CLs_obs, CLs_exp = pyhf.infer.hypotest(
    1.0,  # null hypothesis
    [52.5, 65.0] + model.config.auxdata,
    model,
    test_stat="qtilde",
    return_expected_set=True,
)
print(f"      Observed CLs: {CLs_obs:.4f}")
for expected_value, n_sigma in zip(CLs_exp, np.arange(-2, 3)):
    print(f"Expected CLs({n_sigma:2d} σ): {expected_value:.4f}")<jupyter_output><empty_output><jupyter_text>Simple Upper LimitTo get upper limits, we just need to run multiple hypothesis tests for a lot of different null hypotheses of BSM with $\mu \in [0, \ldots, 5.0]$ and then find the value of $\mu$ for which the null hypothesis is rejected (a 95% $\text{CL}_\text{s}$).<jupyter_code>poi_values = np.linspace(0.1, 5, 50)
results = [
    pyhf.infer.hypotest(
        poi_value,
        observations,
        model,
        test_stat="qtilde",
        return_expected_set=True,
    )
    for poi_value in poi_values
]<jupyter_output><empty_output><jupyter_text>So this ran pretty fast. We can plot the "standard Brazil band" using the `pyhf.contrib` module (which needs `pyhf[contrib]`):<jupyter_code>fig, ax = plt.subplots()
fig.set_size_inches(10.5, 7)
ax.set_title("Hypothesis Tests")
ax.set_xlabel(r"$\mu$")
ax.set_ylabel(r"$\mathrm{CL}_{s}$")

brazil.plot_results(ax, poi_values, results)<jupyter_output><empty_output><jupyter_text>And then just calculate the upper limit by interpolating:<jupyter_code>observed = np.asarray([h[0] for h in results]).ravel()
expected = np.asarray([h[1][2] for h in results]).ravel()
print(f"Upper limit (obs): μ = {np.interp(0.05, observed[::-1], poi_values[::-1]):.4f}")
print(f"Upper limit (exp): μ = {np.interp(0.05, expected[::-1], poi_values[::-1]):.4f}")<jupyter_output><empty_output><jupyter_text>Alternatively, you can also just use the [`upperlimit` API](https://pyhf.readthedocs.io/en/v0.6.0/_generated/pyhf.infer.intervals.upperlimit.html)<jupyter_code>obs_limit, exp_limits, (scan, results) = pyhf.infer.intervals.upperlimit(
    observations, model, poi_values, level=0.05, return_results=True
)
print(f"Upper limit (obs): μ = {obs_limit:.4f}")
print(f"Upper limit (exp): μ = {exp_limits[2]:.4f}")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11934935981186308, 'fraction_numerical': 0.013718317219754377, 'mean_word_length': 5.352282157676348, 'pattern_counts': {'":': 0, '<': 114, '<?xml version=': 0, '>': 113, 'https://': 7, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '28', 'dir': 'jupyter-structured-clean-dedup', 'id': '88088', 'license': 'MIT', 'n_tokens_mistral': 4945, 'n_tokens_neox': 4692, 'n_words': 1631, 'path': 'book/HelloWorld.ipynb', 'repo_name': 'pyhf/tutorial-Reinterpretation-Forum-2021', 'size': '15308'}