hypotests ######### This submodule provides tools to do statistical inferences such as discovery test and computations of upper limits or confidence intervals. hepstats needs a fitting backend to perform computations such as zfit _. Any fitting library can be used if their API is compatible with hepstats (see api checks _). We give here a simple example of an upper limit calculation of the yield of a Gaussian signal with known mean and sigma over an exponential background. The fitting backend used is the zfit _ package. If you are unfamiliar with zfit you can have a look at the zfit documentation _. First we import what's necessary from zfit, such as the **ExtendedUnbinnedNLL** class as we want to construct an extended unbinned likelihood. **Minuit** is also imported, it is a zfit wrapper of the the minuit minimizer from iminuit _. .. code-block:: pycon >>> import zfit >>> from zfit.loss import ExtendedUnbinnedNLL >>> from zfit.minimize import Minuit >>> import numpy as np Then we construct the data sample which consists of 300 points that are drawn from an exponential distribution with -2 slope, constituting the background, whereas 10 points drawn from a Gaussian distribution of mean 1.2 and width 0.1, is the signal. The fit range is defined between 0.1 and 3.0 meaning that some points of the background distribution are filtered out. The data, which is a numpy array, is then transformed into a zfit **Data** object. .. code-block:: pycon >>> bounds = (0.1, 3.0) >>> obs = zfit.Space('x', limits=bounds) >>> bkg = np.random.exponential(1/2, 300) >>> peak = np.random.normal(1.2, 0.1, 10) >>> data = np.concatenate((bkg, peak)) >>> data = data[(data > bounds[0]) & (data < bounds[1])] >>> data = zfit.Data.from_numpy(obs=obs, array=data) Now we build the model. For the background an exponential pdf with **lambda_**, the slope of the exponential as a free parameter. For the signal a Gaussian pdf is used with mean and width fixed to 1.2 and 0.1 respectively. The background and signal pdfs are extended using the yield parameters **Nbkg** and **Nsig** respectively, which are free. The extended negative log-likelihood is then constructed using the background and signal models summed and the data. .. code-block:: pycon >>> lambda_ = zfit.Parameter("lambda", -2.0, -4.0, -1.0) >>> Nsig = zfit.Parameter("Nsig", 1., -20., 500) >>> Nbkg = zfit.Parameter("Nbkg", 250, 0., 50) >>> signal = zfit.pdf.Gauss(obs=obs, mu=1.2, sigma=0.1).create_extended(Nsig) >>> background = zfit.pdf.Exponential(obs=obs, lambda_=lambda_).create_extended(Nbkg) >>> total = zfit.pdf.SumPDF([signal, background]) >>> nll = ExtendedUnbinnedNLL(model=total, data=data) The background plus signal can then be fitted to the data. .. code-block:: pycon >>> # Instantiate a minuit minimizer >>> minimizer = Minuit() >>> # minimisation of the loss function >>> minimum = minimizer.minimize(loss=nll) >>> minimum.hesse() >>> print(minimum) +---------+-------------+------------------+---------+-------------+ | valid | converged | param at limit | edm | min value | +=========+=============+==================+=========+=============+ | True | True | False | 4.9e-05 | -1077 | +---------+-------------+------------------+---------+-------------+ Parameters +--------+--------+---------------+-----------+ | name | value | hesse | at limit | +========+========+===============+===========+ | Nsig | 4.518 | +/- 5.8 | False | +--------+--------+---------------+-----------+ | Nbkg | 251.6 | +/- 17 | False | +--------+--------+---------------+-----------+ | lambda | -1.93 | +/- 0.14 | False | +--------+--------+---------------+-----------+ So the fitted number of signal candidates is 4.518 +/- 5.8, which is consistent with zero. We can then compute an upper limit on this number which should be approximately equal to 4.5 + 2 * 5.8 ≈ 16. First we import from the :py:mod:~hepstats.hypotests.calculators submodule of :py:mod:~hepstats.hypotests the :py:class:~hepstats.hypotests.calculators.asymptotic_calculator.AsymptoticCalculator which takes as input the loss function and minimizer. >>> from hepstats.hypotests.calculators import AsymptoticCalculator >>> calculator = AsymptoticCalculator(nll, Minuit(), asimov_bins=100) The :py:class:~hepstats.hypotests.parameters.POI and :py:class:~hepstats.hypotests.parameters.POIarray classes are also imported, POI stands for parameter of interest. In our case the POI is **Nsig**. To compute an upper limit you need to explicitly specify the background-only hypothesis (null) and the background plus signal hypothesis, in hepstats this done using :py:class:~hepstats.hypotests.parameters.POI/ :py:class:~hepstats.hypotests.parameters.POIarray: .. code-block:: pycon >>> from hepstats.hypotests.parameters import POI, POIarray >>> >>> # background only >>> poialt = POI(Nsig, 0) >>> # background + signal >>> poinull = POIarray(Nsig, np.linspace(0.0, 25, 20)) A :py:class:~hepstats.hypotests.parameters.POI takes as input the parameter **Nsig** and a single value for a given hypothesis, for **poialt** it's 0 because this is the background only hypothesis. Similarly :py:class:~hepstats.hypotests.parameters.POIarray takes as input the parameter **Nsig** and an array of values to scan for **Nsig**, from 0 to 25. A range is needed because the **calculator** instance will compute a *p-value* for each value in **poinull**, the upper limit for a given confidence level :math:\alpha is defined as the value of **Nsig** for which the *p-value* is equal to :math:1 - \alpha. We can now create an :py:class:~hepstats.hypotests.core.upperlimit.UpperLimit instance which takes as input the **calculator**, **poinull** and **poialt**. The :py:class:~hepstats.hypotests.core.upperlimit.UpperLimit instance will ask the **calculator** to compute the *p-values* for each value in **poinull** and eventually find the value of the upper limit on **Nsig** (if the upper limit is in the range of the **poinull** values). Below is an example on how to compute a CLs upper limit at 95 % confidence level. .. code-block:: pycon >>> from hepstats.hypotests import UpperLimit >>> ul = UpperLimit(calculator, poinull, poialt) >>> ul.upperlimit(alpha=0.05, CLs=True) Observed upper limit: Nsig = 15.725784747406346 Expected upper limit: Nsig = 11.927442041887158 Expected upper limit +1 sigma: Nsig = 16.596396280677116 Expected upper limit -1 sigma: Nsig = 8.592750403611896 Expected upper limit +2 sigma: Nsig = 22.24864429383046 Expected upper limit -2 sigma: Nsig = 6.400549971360598 In the result you obtain the observed and expected limits. The observed limit is the limit based on the observation of 4.518 +/- 5.8 signal candidates in data. The expected limit is the limit under the background only hypothesis. A graphical representation on how the upper limit is computed in shown in the following figure. .. image:: https://raw.githubusercontent.com/scikit-hep/hepstats/master/notebooks/hypotests/asy_ul.png