# 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 observations

This is probably the simplest model one could use in HistFactory... so let's make it.


```python
import json
import numpy as np
import matplotlib.pyplot as plt

import pyhf
from pyhf.contrib.viz import brazil  # not imported by default!
```


```python
model = pyhf.simplemodels.hepdata_like(
    signal_data=[5.0, 10.0], bkg_data=[50.0, 60.0], bkg_uncerts=[5.0, 12.0]
)
model
```

What did we just make? This returns a [`pyhf.Model`](https://pyhf.readthedocs.io/en/v0.6.0/_generated/pyhf.pdf.Model.html#pyhf.pdf.Model) object. Let's check out the specification.


```python
print(json.dumps(model.spec, indent=2))
```

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.


```python
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}")
```

## 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 Poisson
2. the auxiliary data is fully determined by the shapesys 'data' and the 'background' data

So 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.html#id24) in the introduction of our pyhf documentation to use as reference.



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:


```python
(np.array([5.0, 12.0]) / np.array([50.0, 60.0])) ** -2
```


```python
model.config.auxdata
```

Whew. Ok, I guess that makes sense.

## Data and Parameters

What about the actual data for the model then?


```python
model.expected_data()
```

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:


```python
model.expected_data([1.0, 1.0, 1.0])
```

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.


```python
model.expected_data([1.0, 1.0, 1.0], include_auxdata=False)
```

Similarly, there are also methods separately for the actual data and the auxdata, both which take in all the parameters:


```python
model.expected_actualdata([1.0, 1.0, 1.0])
```


```python
model.expected_auxdata([1.0, 1.0, 1.0])
```

### 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:


```python
model.config.parameters
```


```python
model.config.modifiers
```

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


```python
model.config.par_order
```

### 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:


```python
model.config.channel_nbins
```

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`:


```python
model.config.param_set("uncorr_bkguncrt")
```

which we can see is constrained by a Poisson and


```python
model.config.param_set("uncorr_bkguncrt").n_parameters
```

### 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...


```python
model.config.suggested_init()
```

... and their bounds...


```python
model.config.suggested_bounds()
```

... and if the parameter should be held constant in a fit or not ...


```python
model.config.suggested_fixed()
```

So armed with this knowledge, we could have also done


```python
init_pars = model.config.suggested_init()
model.expected_actualdata(init_pars)
```

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:


```python
model.config.poi_index
```

and then just change the parameter


```python
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)
```

## Simple Inference

The 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.


```python
observations = [52.5, 65.0] + model.config.auxdata  # this is a common pattern!

model.logpdf(pars=bkg_pars, data=observations)
```


```python
model.logpdf(pars=init_pars, data=observations)
```

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.html#inference) 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.


```python
pyhf.infer.mle.fit(data=observations, pdf=model)
```

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 Testing

Great, 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 fit

So 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 fit

So let's run a hypothesis test for

* null hypothesis ($\mu = 1$) - "SUSY is real"
* alternate hypothesis ($\mu = 0$) - "Standard Model explains it all"


```python
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}")
```

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).


```python
model.config.suggested_bounds()[model.config.poi_index]
```

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$.


```python
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}")
```

## Simple Upper Limit

To 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}$).


```python
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
]
```

So this ran pretty fast. We can plot the "standard Brazil band" using the `pyhf.contrib` module (which needs `pyhf[contrib]`):


```python
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)
```

And then just calculate the upper limit by interpolating:


```python
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}")
```

Alternatively, you can also just use the [`upperlimit` API](https://pyhf.readthedocs.io/en/v0.6.0/_generated/pyhf.infer.intervals.upperlimit.html)


```python
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}")
```

4171283443dcf359aa55c04a18a00683d08098f9
ipynb
Jupyter Notebook
book/HelloWorld.ipynb
pyhf/tutorial-Reinterpretation-Forum-2021
b1136f98a8257cd20207d2da5a14a848d3a97437
MIT
2021-02-18T12:46:49.000Z
2021-02-18T12:46:49.000Z
book/HelloWorld.ipynb
pyhf/tutorial-Reinterpretation-Forum-2021
b1136f98a8257cd20207d2da5a14a848d3a97437
MIT
2021-02-19T00:24:27.000Z
2021-02-19T00:24:27.000Z
book/HelloWorld.ipynb
pyhf/tutorial-Reinterpretation-Forum-2021
b1136f98a8257cd20207d2da5a14a848d3a97437
MIT
2021-02-23T20:17:24.000Z
2021-02-23T20:17:24.000Z
Qwen/Qwen-72B
1. YES
2. YES