Introduction

Measurements in High Energy Physics (HEP) rely on determining the compatibility of observed collision events with theoretical predictions. The relationship between them is often formalised in a statistical model \(f(\bm{x}|\fullset)\) describing the probability of data \(\bm{x}\) given model parameters \(\fullset\). Given observed data, the likelihood \(\mathcal{L}(\fullset)\) then serves as the basis to test hypotheses on the parameters \(\fullset\). For measurements based on binned data (histograms), the \(\HiFa{}\) family of statistical models has been widely used in both Standard Model measurements [intro-4] as well as searches for new physics [intro-5]. In this package, a declarative, plain-text format for describing \(\HiFa{}\)-based likelihoods is presented that is targeted for reinterpretation and long-term preservation in analysis data repositories such as HEPData [intro-3].

HistFactory

Statistical models described using \(\HiFa{}\) [intro-2] center around the simultaneous measurement of disjoint binned distributions (channels) observed as event counts \(\channelcounts\). For each channel, the overall expected event rate 1 is the sum over a number of physics processes (samples). The sample rates may be subject to parametrised variations, both to express the effect of free parameters \(\freeset\) 2 and to account for systematic uncertainties as a function of constrained parameters \(\constrset\). The degree to which the latter can cause a deviation of the expected event rates from the nominal rates is limited by constraint terms. In a frequentist framework these constraint terms can be viewed as auxiliary measurements with additional global observable data \(\auxdata\), which paired with the channel data \(\channelcounts\) completes the observation \(\bm{x} = (\channelcounts,\auxdata)\). In addition to the partition of the full parameter set into free and constrained parameters \(\fullset = (\freeset,\constrset)\), a separate partition \(\fullset = (\poiset,\nuisset)\) will be useful in the context of hypothesis testing, where a subset of the parameters are declared parameters of interest \(\poiset\) and the remaining ones as nuisance parameters \(\nuisset\).

(1)\[f(\bm{x}|\fullset) = f(\bm{x}|\overbrace{\freeset}^{\llap{\text{free}}},\underbrace{\constrset}_{\llap{\text{constrained}}}) = f(\bm{x}|\overbrace{\poiset}^{\rlap{\text{parameters of interest}}},\underbrace{\nuisset}_{\rlap{\text{nuisance parameters}}})\]

Thus, the overall structure of a \(\HiFa{}\) probability model is a product of the analysis-specific model term describing the measurements of the channels and the analysis-independent set of constraint terms:

(2)\[\begin{split}f(\channelcounts, \auxdata \,|\,\freeset,\constrset) = \underbrace{\color{blue}{\prod_{c\in\mathrm{\,channels}} \prod_{b \in \mathrm{\,bins}_c}\textrm{Pois}\left(n_{cb} \,\middle|\, \nu_{cb}\left(\freeset,\constrset\right)\right)}}_{\substack{\text{Simultaneous measurement}\\% \text{of multiple channels}}} \underbrace{\color{red}{\prod_{\singleconstr \in \constrset} c_{\singleconstr}(a_{\singleconstr} |\, \singleconstr)}}_{\substack{\text{constraint terms}\\% \text{for }\unicode{x201C}\text{auxiliary measurements}\unicode{x201D}}},\end{split}\]

where within a certain integrated luminosity we observe \(n_{cb}\) events given the expected rate of events \(\nu_{cb}(\freeset,\constrset)\) as a function of unconstrained parameters \(\freeset\) and constrained parameters \(\constrset\). The latter has corresponding one-dimensional constraint terms \(c_\singleconstr(a_\singleconstr|\,\singleconstr)\) with auxiliary data \(a_\singleconstr\) constraining the parameter \(\singleconstr\). The event rates \(\nu_{cb}\) are defined as

(3)\[\nu_{cb}\left(\fullset\right) = \sum_{s\in\mathrm{\,samples}} \nu_{scb}\left(\freeset,\constrset\right) = \sum_{s\in\mathrm{\,samples}}\underbrace{\left(\prod_{\kappa\in\,\bm{\kappa}} \kappa_{scb}\left(\freeset,\constrset\right)\right)}_{\text{multiplicative modifiers}}\, \Bigg(\nu_{scb}^0\left(\freeset, \constrset\right) + \underbrace{\sum_{\Delta\in\bm{\Delta}} \Delta_{scb}\left(\freeset,\constrset\right)}_{\text{additive modifiers}}\Bigg)\,.\]

The total rates are the sum over sample rates \(\nu_{csb}\), each determined from a nominal rate \(\nu_{scb}^0\) and a set of multiplicative and additive denoted rate modifiers \(\bm{\kappa}(\fullset)\) and \(\bm{\Delta}(\fullset)\). These modifiers are functions of (usually a single) model parameters. Starting from constant nominal rates, one can derive the per-bin event rate modification by iterating over all sample rate modifications as shown in (3).

As summarised in Modifiers and Constraints, rate modifications are defined in \(\HiFa{}\) for bin \(b\), sample \(s\), channel \(c\). Each modifier is represented by a parameter \(\phi \in \{\gamma, \alpha, \lambda, \mu\}\). By convention bin-wise parameters are denoted with \(\gamma\) and interpolation parameters with \(\alpha\). The luminosity \(\lambda\) and scale factors \(\mu\) affect all bins equally. For constrained modifiers, the implied constraint term is given as well as the necessary input data required to construct it. \(\sigma_b\) corresponds to the relative uncertainty of the event rate, whereas \(\delta_b\) is the event rate uncertainty of the sample relative to the total event rate \(\nu_b = \sum_s \nu^0_{sb}\).

Modifiers implementing uncertainties are paired with a corresponding default constraint term on the parameter limiting the rate modification. The available modifiers may affect only the total number of expected events of a sample within a given channel, i.e. only change its normalisation, while holding the distribution of events across the bins of a channel, i.e. its “shape”, invariant. Alternatively, modifiers may change the sample shapes. Here \(\HiFa{}\) supports correlated an uncorrelated bin-by-bin shape modifications. In the former, a single nuisance parameter affects the expected sample rates within the bins of a given channel, while the latter introduces one nuisance parameter for each bin, each with their own constraint term. For the correlated shape and normalisation uncertainties, \(\HiFa{}\) makes use of interpolating functions, \(f_p\) and \(g_p\), constructed from a small number of evaluations of the expected rate at fixed values of the parameter \(\alpha\) 3. For the remaining modifiers, the parameter directly affects the rate.

Modifiers and Constraints

Description	Modification	Constraint Term \(c_\singleconstr\)	Input
Uncorrelated Shape	\(\kappa_{scb}(\gamma_b) = \gamma_b\)	\(\prod_b \mathrm{Pois}\left(r_b = \sigma_b^{-2}\middle|\,\rho_b = \sigma_b^{-2}\gamma_b\right)\)	\(\sigma_{b}\)
Correlated Shape	\(\Delta_{scb}(\alpha) = f_p\left(\alpha\middle|\,\Delta_{scb,\alpha=-1},\Delta_{scb,\alpha = 1}\right)\)	\(\displaystyle\mathrm{Gaus}\left(a = 0\middle|\,\alpha,\sigma = 1\right)\)	\(\Delta_{scb,\alpha=\pm1}\)
Normalisation Unc.	\(\kappa_{scb}(\alpha) = g_p\left(\alpha\middle|\,\kappa_{scb,\alpha=-1},\kappa_{scb,\alpha=1}\right)\)	\(\displaystyle\mathrm{Gaus}\left(a = 0\middle|\,\alpha,\sigma = 1\right)\)	\(\kappa_{scb,\alpha=\pm1}\)
MC Stat. Uncertainty	\(\kappa_{scb}(\gamma_b) = \gamma_b\)	\(\prod_b \mathrm{Gaus}\left(a_{\gamma_b} = 1\middle|\,\gamma_b,\delta_b\right)\)	\(\delta_b^2 = \sum_s\delta^2_{sb}\)
Luminosity	\(\kappa_{scb}(\lambda) = \lambda\)	\(\displaystyle\mathrm{Gaus}\left(l = \lambda_0\middle|\,\lambda,\sigma_\lambda\right)\)	\(\lambda_0,\sigma_\lambda\)
Normalisation	\(\kappa_{scb}(\mu_b) = \mu_b\)
Data-driven Shape	\(\kappa_{scb}(\gamma_b) = \gamma_b\)

Given the likelihood \(\mathcal{L}(\fullset)\), constructed from observed data in all channels and the implied auxiliary data, measurements in the form of point and interval estimates can be defined. The majority of the parameters are nuisance parameters — parameters that are not the main target of the measurement but are necessary to correctly model the data. A small subset of the unconstrained parameters may be declared as parameters of interest for which measurements hypothesis tests are performed, e.g. profile likelihood methods [intro-1]. The Symbol Notation table provides a summary of all the notation introduced in this documentation.

Symbol Notation

Symbol	Name
\(f(\bm{x} | \fullset)\)	model
\(\mathcal{L}(\fullset)\)	likelihood
\(\bm{x} = \{\channelcounts, \auxdata\}\)	full dataset (including auxiliary data)
\(\channelcounts\)	channel data (or event counts)
\(\auxdata\)	auxiliary data
\(\nu(\fullset)\)	calculated event rates
\(\fullset = \{\freeset, \constrset\} = \{\poiset, \nuisset\}\)	all parameters
\(\freeset\)	free parameters
\(\constrset\)	constrained parameters
\(\poiset\)	parameters of interest
\(\nuisset\)	nuisance parameters
\(\bm{\kappa}(\fullset)\)	multiplicative rate modifier
\(\bm{\Delta}(\fullset)\)	additive rate modifier
\(c_\singleconstr(a_\singleconstr | \singleconstr)\)	constraint term for constrained parameter \(\singleconstr\)
\(\sigma_\singleconstr\)	relative uncertainty in the constrained parameter

Declarative Formats

While flexible enough to describe a wide range of LHC measurements, the design of the \(\HiFa{}\) specification is sufficiently simple to admit a declarative format that fully encodes the statistical model of the analysis. This format defines the channels, all associated samples, their parameterised rate modifiers and implied constraint terms as well as the measurements. Additionally, the format represents the mathematical model, leaving the implementation of the likelihood minimisation to be analysis-dependent and/or language-dependent. Originally XML was chosen as a specification language to define the structure of the model while introducing a dependence on \(\Root{}\) to encode the nominal rates and required input data of the constraint terms [intro-2]. Using this specification, a model can be constructed and evaluated within the \(\RooFit{}\) framework.

This package introduces an updated form of the specification based on the ubiquitous plain-text JSON format and its schema-language JSON Schema. Described in more detail in Likelihood Specification, this schema fully specifies both structure and necessary constrained data in a single document and thus is implementation independent.

Additional Material

Footnotes

1

Here rate refers to the number of events expected to be observed within a given data-taking interval defined through its integrated luminosity. It often appears as the input parameter to the Poisson distribution, hence the name “rate”.

2

These free parameters frequently include the of a given process, i.e. its cross-section normalised to a particular reference cross-section such as that expected from the Standard Model or a given BSM scenario.

3

This is usually constructed from the nominal rate and measurements of the event rate at \(\alpha=\pm1\), where the value of the modifier at \(\alpha=\pm1\) must be provided and the value at \(\alpha=0\) corresponds to the corresponding identity operation of the modifier, i.e. \(f_{p}(\alpha=0) = 0\) and \(g_{p}(\alpha = 0)=1\) for additive and multiplicative modifiers respectively. See Section 4.1 in [intro-2].

Bibliography

intro-1

Glen Cowan, Kyle Cranmer, Eilam Gross, and Ofer Vitells. Asymptotic formulae for likelihood-based tests of new physics. Eur. Phys. J. C, 71:1554, 2011. arXiv:1007.1727, doi:10.1140/epjc/s10052-011-1554-0.

intro-2(1,2,3)

Kyle Cranmer, George Lewis, Lorenzo Moneta, Akira Shibata, and Wouter Verkerke. HistFactory: A tool for creating statistical models for use with RooFit and RooStats. Technical Report CERN-OPEN-2012-016, New York U., New York, Jan 2012. URL: https://cds.cern.ch/record/1456844.

intro-3

Eamonn Maguire, Lukas Heinrich, and Graeme Watt. HEPData: a repository for high energy physics data. J. Phys. Conf. Ser., 898(10):102006, 2017. arXiv:1704.05473, doi:10.1088/1742-6596/898/10/102006.

intro-4

ATLAS Collaboration. Measurements of Higgs boson production and couplings in diboson final states with the ATLAS detector at the LHC. Phys. Lett. B, 726:88, 2013. arXiv:1307.1427, doi:10.1016/j.physletb.2014.05.011.

intro-5

ATLAS Collaboration. Search for supersymmetry in final states with missing transverse momentum and multiple \(b\)-jets in proton–proton collisions at \(\sqrt s = 13\) \(\TeV \) with the ATLAS detector. ATLAS-CONF-2018-041, 2018. URL: https://cds.cern.ch/record/2632347.

Likelihood Specification

The structure of the JSON specification of models follows closely the original XML-based specification [likelihood-2].

Workspace

{
    "$schema": "http://json-schema.org/draft-06/schema#",
    "$id": "https://scikit-hep.org/pyhf/schemas/1.0.0/workspace.json",
    "$ref": "defs.json#/definitions/workspace"
}

The overall document in the above code snippet describes a workspace, which includes

• channels: The channels in the model, which include a description of the samples within each channel and their possible parametrised modifiers.

• measurements: A set of measurements, which define among others the parameters of interest for a given statistical analysis objective.

• observations: The observed data, with which a likelihood can be constructed from the model.

A workspace consists of the channels, one set of observed data, but can include multiple measurements. If provided a JSON file, one can quickly check that it conforms to the provided workspace specification as follows:

import json, requests, jsonschema

workspace = json.load(open("/path/to/analysis_workspace.json"))
# if no exception is raised, it found and parsed the schema
schema = requests.get("https://scikit-hep.org/pyhf/schemas/1.0.0/workspace.json").json()
# If no exception is raised by validate(), the instance is valid.
jsonschema.validate(instance=workspace, schema=schema)

Channel

A channel is defined by a channel name and a list of samples [likelihood-1].

{
    "channel": {
        "type": "object",
        "properties": {
            "name": { "type": "string" },
            "samples": { "type": "array", "items": {"$ref": "#/definitions/sample"}, "minItems": 1 }
        },
        "required": ["name", "samples"],
        "additionalProperties": false
    },
}

The Channel specification consists of a list of channel descriptions. Each channel, an analysis region encompassing one or more measurement bins, consists of a name field and a samples field (see Channel), which holds a list of sample definitions (see Sample). Each sample definition in turn has a name field, a data field for the nominal event rates for all bins in the channel, and a modifiers field of the list of modifiers for the sample.

Sample

A sample is defined by a sample name, the sample event rate, and a list of modifiers [likelihood-1].

{
    "sample": {
        "type": "object",
        "properties": {
            "name": { "type": "string" },
            "data": { "type": "array", "items": {"type": "number"}, "minItems": 1 },
            "modifiers": {
                "type": "array",
                "items": {
                    "anyOf": [
                        { "$ref": "#/definitions/modifier/histosys" },
                        { "$ref": "#/definitions/modifier/lumi" },
                        { "$ref": "#/definitions/modifier/normfactor" },
                        { "$ref": "#/definitions/modifier/normsys" },
                        { "$ref": "#/definitions/modifier/shapefactor" },
                        { "$ref": "#/definitions/modifier/shapesys" },
                        { "$ref": "#/definitions/modifier/staterror" }
                    ]
                }
            }
        },
        "required": ["name", "data", "modifiers"],
        "additionalProperties": false
    },
}

Modifiers

The modifiers that are applicable for a given sample are encoded as a list of JSON objects with three fields. A name field, a type field denoting the class of the modifier, and a data field which provides the necessary input data as denoted in Modifiers and Constraints.

Based on the declared modifiers, the set of parameters and their constraint terms are derived implicitly as each type of modifier unambiguously defines the constraint terms it requires. Correlated shape modifiers and normalisation uncertainties have compatible constraint terms and thus modifiers can be declared that share parameters by re-using a name 1 for multiple modifiers. That is, a variation of a single parameter causes a shift within sample rates due to both shape and normalisation variations.

We review the structure of each modifier type below.

Uncorrelated Shape (shapesys)

To construct the constraint term, the relative uncertainties \(\sigma_b\) are necessary for each bin. Therefore, we record the absolute uncertainty as an array of floats, which combined with the nominal sample data yield the desired \(\sigma_b\). An example is shown below:

{ "name": "mod_name", "type": "shapesys", "data": [1.0, 1.5, 2.0] }

An example of an uncorrelated shape modifier with three absolute uncertainty terms for a 3-bin channel.

Warning

Nuisance parameters will not be allocated for any bins where either

• the samples nominal expected rate is zero, or

• the absolute uncertainty is zero.

These values are, in the context of uncorrelated shape uncertainties, unphysical. If this situation occurs, one needs to go back and understand the inputs as this is undefined behavior in HistFactory.

The previous example will allocate three nuisance parameters for mod_name. The following example will allocate only two nuisance parameters for a 3-bin channel:

{ "name": "mod_name", "type": "shapesys", "data": [1.0, 0.0, 2.0] }

Correlated Shape (histosys)

This modifier represents the same source of uncertainty which has a different effect on the various sample shapes, hence a correlated shape. To implement an interpolation between sample distribution shapes, the distributions with a “downward variation” (“lo”) associated with \(\alpha=-1\) and an “upward variation” (“hi”) associated with \(\alpha=+1\) are provided as arrays of floats. An example is shown below:

{ "name": "mod_name", "type": "histosys", "data": {"hi_data": [20,15], "lo_data": [10, 10]} }

An example of a correlated shape modifier with absolute shape variations for a 2-bin channel.

Normalisation Uncertainty (normsys)

The normalisation uncertainty modifies the sample rate by a overall factor \(\kappa(\alpha)\) constructed as the interpolation between downward (“lo”) and upward (“hi”) as well as the nominal setting, i.e. \(\kappa(-1) = \kappa_{\alpha=-1}\), \(\kappa(0) = 1\) and \(\kappa(+1) = \kappa_{\alpha=+1}\). In the modifier definition we record \(\kappa_{\alpha=+1}\) and \(\kappa_{\alpha=-1}\) as floats. An example is shown below:

{ "name": "mod_name", "type": "normsys", "data": {"hi": 1.1, "lo": 0.9} }

An example of a normalisation uncertainty modifier with scale factors recorded for the up/down variations of an \(n\)-bin channel.

MC Statistical Uncertainty (staterror)

As the sample counts are often derived from Monte Carlo (MC) datasets, they necessarily carry an uncertainty due to the finite sample size of the datasets. As explained in detail in [likelihood-2], adding uncertainties for each sample would yield a very large number of nuisance parameters with limited utility. Therefore a set of bin-wise scale factors \(\gamma_b\) is introduced to model the overall uncertainty in the bin due to MC statistics. The constrained term is constructed as a set of Gaussian constraints with a central value equal to unity for each bin in the channel. The scales \(\sigma_b\) of the constraint are computed from the individual uncertainties of samples defined within the channel relative to the total event rate of all samples: \(\delta_{csb} = \sigma_{csb}/\sum_s \nu^0_{scb}\). As not all samples are within a channel are estimated from MC simulations, only the samples with a declared statistical uncertainty modifier enter the sum. An example is shown below:

{ "name": "mod_name", "type": "staterror", "data": [0.1] }

An example of a statistical uncertainty modifier.

Luminosity (lumi)

Sample rates derived from theory calculations, as opposed to data-driven estimates, are scaled to the integrated luminosity corresponding to the observed data. As the luminosity measurement is itself subject to an uncertainty, it must be reflected in the rate estimates of such samples. As this modifier is of global nature, no additional per-sample information is required and thus the data field is nulled. This uncertainty is relevant, in particular, when the parameter of interest is a signal cross-section. The luminosity uncertainty \(\sigma_\lambda\) is provided as part of the parameter configuration included in the measurement specification discussed in Measurements. An example is shown below:

{ "name": "mod_name", "type": "lumi", "data": null }

An example of a luminosity modifier.

Unconstrained Normalisation (normfactor)

The unconstrained normalisation modifier scales the event rates of a sample by a free parameter \(\mu\). Common use cases are the signal rate of a possible BSM signal or simultaneous in-situ measurements of background samples. Such parameters are frequently the parameters of interest of a given measurement. No additional per-sample data is required. An example is shown below:

{ "name": "mod_name", "type": "normfactor", "data": null }

An example of a normalisation modifier.

Data-driven Shape (shapefactor)

In order to support data-driven estimation of sample rates (e.g. for multijet backgrounds), the data-driven shape modifier adds free, bin-wise multiplicative parameters. Similarly to the normalisation factors, no additional data is required as no constraint is defined. An example is shown below:

{ "name": "mod_name", "type": "shapefactor", "data": null }

An example of an uncorrelated shape modifier.

Data

The data provided by the analysis are the observed data for each channel (or region). This data is provided as a mapping from channel name to an array of floats, which provide the observed rates in each bin of the channel. The auxiliary data is not included as it is an input to the likelihood that does not need to be archived and can be determined automatically from the specification. An example is shown below:

{ "chan_name_one": [10, 20], "chan_name_two": [4, 0]}

An example of channel data.

Measurements

Given the data and the model definitions, a measurement can be defined. In the current schema, the measurements defines the name of the parameter of interest as well as parameter set configurations. 2 Here, the remaining information not covered through the channel definition is provided, e.g. for the luminosity parameter. For all modifiers, the default settings can be overridden where possible:

• inits: Initial value of the parameter.

• bounds: Interval bounds of the parameter.

• auxdata: Auxiliary data for the associated constraint term.

• sigmas: Associated uncertainty of the parameter.

An example is shown below:

{
    "name": "MyMeasurement",
    "config": {
        "poi": "SignalCrossSection", "parameters": [
            { "name":"lumi", "auxdata":[1.0],"sigmas":[0.017], "bounds":[[0.915,1.085]],"inits":[1.0] },
            { "name":"mu_ttbar", "bounds":[[0, 5]] },
            { "name":"rw_1CR", "fixed":true }
        ]
    }
}

An example of a measurement. This measurement, which scans over the parameter of interest SignalCrossSection, is setting configurations for the luminosity modifier, changing the default bounds for the normfactor modifier named mu_ttbar, and specifying that the modifier rw_1CR is held constant (fixed).

Observations

This is what we evaluate the hypothesis testing against, to determine the compatibility of signal+background hypothesis to the background-only hypothesis. This is specified as a list of objects, with each object structured as

• name: the channel for which the observations are recorded

• data: the bin-by-bin observations for the named channel

An example is shown below:

{
    "name": "channel1",
    "data": [110.0, 120.0]
}

An example of an observation. This observation recorded for a 2-bin channel channel1, has values 110.0 and 120.0.

Toy Example

{
    "channels": [
        { "name": "singlechannel",
          "samples": [
            { "name": "signal",
              "data": [5.0, 10.0],
              "modifiers": [ { "name": "mu", "type": "normfactor", "data": null} ]
            },
            { "name": "background",
              "data": [50.0, 60.0],
              "modifiers": [ {"name": "uncorr_bkguncrt", "type": "shapesys", "data": [5.0, 12.0]} ]
            }
          ]
        }
    ],
    "observations": [
        { "name": "singlechannel", "data": [50.0, 60.0] }
    ],
    "measurements": [
        { "name": "Measurement", "config": {"poi": "mu", "parameters": []} }
    ],
    "version": "1.0.0"
}

In the above example, we demonstrate a simple measurement of a single two-bin channel with two samples: a signal sample and a background sample. The signal sample has an unconstrained normalisation factor \(\mu\), 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.

Additional Material

Footnotes

1

The name of a modifier specifies the parameter set it is controlled by. Modifiers with the same name share parameter sets.

2

In this context a parameter set corresponds to a named lower-dimensional subspace of the full parameters \(\fullset\). In many cases these are one-dimensional subspaces, e.g. a specific interpolation parameter \(\alpha\) or the luminosity parameter \(\lambda\). For multi-bin channels, however, e.g. all bin-wise nuisance parameters of the uncorrelated shape modifiers are grouped under a single name. Therefore in general a parameter set definition provides arrays of initial values, bounds, etc.

Bibliography

likelihood-1(1,2)

Histfactory definitions schema. Accessed: 2019-06-20. URL: https://scikit-hep.org/pyhf/schemas/1.0.0/defs.json.

likelihood-2(1,2)

Kyle Cranmer, George Lewis, Lorenzo Moneta, Akira Shibata, and Wouter Verkerke. HistFactory: A tool for creating statistical models for use with RooFit and RooStats. Technical Report CERN-OPEN-2012-016, New York U., New York, Jan 2012. URL: https://cds.cern.ch/record/1456844.

Fundamentals

Notebooks:

%pylab inline
from ipywidgets import interact
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

Populating the interactive namespace from numpy and matplotlib

Piecewise Linear Interpolation

References: https://cds.cern.ch/record/1456844/files/CERN-OPEN-2012-016.pdf

We wish to understand interpolation using the piecewise linear function. This is interpcode=0 in the above reference. This function is defined as (nb: vector denotes bold)

\[\eta_s (\vec{\alpha}) = \sigma_{sb}^0(\vec{\alpha}) + \underbrace{\sum_{p \in \text{Syst}} I_\text{lin.} (\alpha_p; \sigma_{sb}^0, \sigma_{psb}^+, \sigma_{psb}^-)}_\text{deltas to calculate}\]

with

\[\begin{split}I_\text{lin.}(\alpha; I^0, I^+, I^-) = \begin{cases} \alpha(I^+ - I^0) \qquad \alpha \geq 0\\ \alpha(I^0 - I^-) \qquad \alpha < 0 \end{cases}\end{split}\]

In this notebook, we’ll demonstrate the technical implementation of these interplations starting from simple dimensionality and increasing the dimensions as we go along. In all situations, we’ll consider a single systematic that we wish to interpolate, such as Jet Energy Scale (JES).

Let’s define the interpolate function. This function will produce the deltas we would like to calculate and sum with the nominal measurement to determine the interpolated measurements value.

def interpolate_deltas(down, nom, up, alpha):
    delta_up = up - nom
    delta_down = nom - down
    if alpha > 0:
        return delta_up * alpha
    else:
        return delta_down * alpha

Why are we calculating deltas? This is some additional foresight that you, the reader, may not have yet. Multiple interpolation schemes exist but they all rely on calculating the change with respect to the nominal measurement (the delta).

Case 1: The Single-binned Histogram

Let’s first start with considering evaluating the total number of events after applying JES corrections. This is the single-bin case. Code that runs through event selection will vary the JES parameter and provide three histograms, each with a single bin. These three histograms represent the nominal-, up-, and down- variations of the JES nuisance parameter.

When processing, we find that there are 10 events nominally, and when we vary the JES parameter downwards, we only measure 8 events. When varying upwards, we measure 15 events.

down_1 = np.array([8])
nom_1 = np.array([10])
up_1 = np.array([15])

We would like to generate a function \(f(\alpha_\text{JES})\) that linearly interpolates the number of events for us so we can scan the phase-space for calculating PDFs. The interpolate_deltas() function defined above does this for us.

alphas = np.linspace(-1.0, 1.0)
deltas = [interpolate_deltas(down_1, nom_1, up_1, alpha) for alpha in alphas]
deltas[:5]

[array([-2.]),
 array([-1.91836735]),
 array([-1.83673469]),
 array([-1.75510204]),
 array([-1.67346939])]

So now that we’ve generated the deltas from the nominal measurement, we can plot this to see how the linear interpolation works in the single-bin case, where we plot the measured values in black, and the interpolation in dashed, blue.

plt.plot(alphas, [nom_1 + delta for delta in deltas], linestyle='--')
plt.scatter((-1, 0, 1), (down_1, nom_1, up_1), color='k')
plt.xlabel(r'$\alpha_\mathrm{JES}$')
plt.ylabel(r'Events')

Text(0,0.5,'Events')

Here, we can imagine building a 1-dimensional tensor (column-vector) of measurements as a function of \(\alpha_\text{JES}\) with each row in the column vector corresponding to a given \(\alpha_\text{JES}\) value.

Case 2: The Multi-binned Histogram

Now, let’s increase the computational difficulty a little by increasing the dimensionality. Assume instead of a single-bin measurement, we have more measurements! We are good physicists after all. Imagine continuing on the previous example, where we add more bins, perhaps because we got more data. Imagine that this was binned by collection year, where we observed 10 events in the first year, 10.5 the next year, and so on…

down_hist = np.linspace(8, 10, 11)
nom_hist = np.linspace(10, 13, 11)
up_hist = np.linspace(15, 20, 11)

Now, we still need to interpolate. Just like before, we have varied JES upwards and downwards to determine the corresponding histograms of variations. In order to interpolate, we need to interpolate by bin for each bin in the three histograms we have here (or three measurements if you prefer).

Let’s go ahead and plot these histograms as a function of the bin index with black as the nominal measurements, red and blue as the down and up variations respectively. The black points are the measurements we have, and for each bin, we would like to interpolate to get an interpolated histogram that represents the measurement as a function of \(\alpha_\text{JES}\).

def plot_measurements(down_hist, nom_hist, up_hist):
    bincenters = np.arange(len(nom_hist))
    for i, h in enumerate(zip(up_hist, nom_hist, down_hist)):
        plt.scatter([i] * len(h), h, color='k', alpha=0.5)

    for c, h in zip(['r', 'k', 'b'], [down_hist, nom_hist, up_hist]):
        plt.plot(bincenters, h, color=c, linestyle='-', alpha=0.5)

    plt.xlabel('Bin index in histogram')
    plt.ylabel('Events')


plot_measurements(down_hist, nom_hist, up_hist)

What does this look like if we evaluate at a single \(\alpha_\text{JES} = 0.5\)? We’ll write a function that interpolates and then plots the interpolated values as a function of bin index, in green, dashed.

def plot_interpolated_histogram(alpha, down_hist, nom_hist, up_hist):
    bincenters = np.arange(len(nom_hist))
    interpolated_vals = [
        nominal + interpolate_deltas(down, nominal, up, alpha)
        for down, nominal, up in zip(down_hist, nom_hist, up_hist)
    ]

    plot_measurements(down_hist, nom_hist, up_hist)
    plt.plot(bincenters, interpolated_vals, color='g', linestyle='--')


plot_interpolated_histogram(0.5, down_hist, nom_hist, up_hist)

We can go one step further in visualization and see what it looks like for different \(\alpha_\text{JES}\) using iPyWidget’s interactivity. Change the slider to get an idea of how the interpolation works.

x = interact(
    lambda alpha: plot_interpolated_histogram(alpha, down_hist, nom_hist, up_hist),
    alpha=(-1, 1, 0.1),
)

The magic in plot_interpolated_histogram() happens to be that for a given \(\alpha_\text{JES}\), we iterate over all measurements bin-by-bin to calculate the interpolated value

[nominal + interpolate_deltas(down, nominal, up, alpha) for down, nominal, up in zip(...hists...)]

So you can imagine that we’re building up a 2-dimensional tensor with each row corresponding to a different \(\alpha_\text{JES}\) and each column corresponding to the bin index of the histograms (or measurements). Let’s go ahead and build a 3-dimensional representation of our understanding so far!

def interpolate_alpha_range(alphas, down_hist, nom_hist, up_hist):
    at_alphas = []
    for alpha in alphas:
        interpolated_hist_at_alpha = [
            nominal + interpolate_deltas(down, nominal, up, alpha)
            for down, nominal, up in zip(down_hist, nom_hist, up_hist)
        ]
        at_alphas.append(interpolated_hist_at_alpha)
    return np.array(at_alphas)

And then with this, we are interpolating over all histograms bin-by-bin and producing a 2-dimensional tensor with each row corresponding to a specific value of \(\alpha_\text{JES}\).

alphas = np.linspace(-1, 1, 11)

interpolated_vals_at_alphas = interpolate_alpha_range(
    alphas, down_hist, nom_hist, up_hist
)

print(interpolated_vals_at_alphas[alphas == -1])
print(interpolated_vals_at_alphas[alphas == 0])
print(interpolated_vals_at_alphas[alphas == 1])

[[ 8.   8.2  8.4  8.6  8.8  9.   9.2  9.4  9.6  9.8 10. ]]
[[10.  10.3 10.6 10.9 11.2 11.5 11.8 12.1 12.4 12.7 13. ]]
[[15.  15.5 16.  16.5 17.  17.5 18.  18.5 19.  19.5 20. ]]

We have a way to generate the 2-dimensional tensor. Let’s go ahead and add in all dimensions. Additionally, we’ll add in some extra code to show the projection of the 2-d plots that we made earlier to help understand the 3-d plot a bit better. Like before, let’s plot specifically colored lines for \(\alpha_\text{JES}=0.5\) as well as provide an interactive session.

def plot_wire(alpha):
    alphas = np.linspace(-1, 1, 51)
    at_alphas = interpolate_alpha_range(alphas, down_hist, nom_hist, up_hist)
    bincenters = np.arange(len(nom_hist))
    x, y = np.meshgrid(bincenters, alphas)
    z = np.asarray(at_alphas)
    bottom = np.zeros_like(x)
    fig = plt.figure(figsize=(10, 10))
    ax1 = fig.add_subplot(111, projection='3d')
    ax1.plot_wireframe(x, y, z, alpha=0.3)

    x, y = np.meshgrid(bincenters, [alpha])
    z = interpolate_alpha_range([alpha], down_hist, nom_hist, up_hist)

    ax1.plot_wireframe(x, y, z, edgecolor='g', linestyle='--')
    ax1.set_xlim(0, 10)
    ax1.set_ylim(-1.0, 1.5)
    ax1.set_zlim(0, 25)
    ax1.view_init(azim=-125)
    ax1.set_xlabel('Bin Index')
    ax1.set_ylabel(r'$\alpha_\mathrm{JES}$')
    ax1.set_zlabel('Events')

    # add in 2D plot goodness

    for c, h, zs in zip(
        ['r', 'k', 'b'], [down_hist, nom_hist, up_hist], [-1.0, 0.0, 1.0]
    ):
        ax1.plot(bincenters, h, color=c, linestyle='-', alpha=0.5, zdir='y', zs=zs)
        ax1.plot(bincenters, h, color=c, linestyle='-', alpha=0.25, zdir='y', zs=1.5)

    ax1.plot(bincenters, z.T, color='g', linestyle='--', zdir='y', zs=alpha)
    ax1.plot(bincenters, z.T, color='g', linestyle='--', alpha=0.5, zdir='y', zs=1.5)

    plt.show()


plot_wire(0.5)

interact(plot_wire, alpha=(-1, 1, 0.1))

<function __main__.plot_wire>

Tensorizing Interpolators

This notebook will introduce some tensor algebra concepts about being able to convert from calculations inside for-loops into a single calculation over the entire tensor. It is assumed that you have some familiarity with what interpolation functions are used for in pyhf.

To get started, we’ll load up some functions we wrote whose job is to generate sets of histograms and alphas that we will compute interpolations for. This allows us to generate random, structured input data that we can use to test the tensorized form of the interpolation function against the original one we wrote. For now, we will consider only the numpy backend for simplicity, but can replace np to pyhf.tensorlib to achieve identical functionality.

The function random_histosets_alphasets_pair will produce a pair (histogramsets, alphasets) of histograms and alphas for those histograms that represents the type of input we wish to interpolate on.

import numpy as np


def random_histosets_alphasets_pair(
    nsysts=150, nhistos_per_syst_upto=300, nalphas=1, nbins_upto=1
):
    def generate_shapes(histogramssets, alphasets):
        h_shape = [len(histogramssets), 0, 0, 0]
        a_shape = (len(alphasets), max(map(len, alphasets)))
        for hs in histogramssets:
            h_shape[1] = max(h_shape[1], len(hs))
            for h in hs:
                h_shape[2] = max(h_shape[2], len(h))
                for sh in h:
                    h_shape[3] = max(h_shape[3], len(sh))
        return tuple(h_shape), a_shape

    def filled_shapes(histogramssets, alphasets):
        # pad our shapes with NaNs
        histos, alphas = generate_shapes(histogramssets, alphasets)
        histos, alphas = np.ones(histos) * np.nan, np.ones(alphas) * np.nan
        for i, syst in enumerate(histogramssets):
            for j, sample in enumerate(syst):
                for k, variation in enumerate(sample):
                    histos[i, j, k, : len(variation)] = variation
        for i, alphaset in enumerate(alphasets):
            alphas[i, : len(alphaset)] = alphaset
        return histos, alphas

    nsyst_histos = np.random.randint(1, 1 + nhistos_per_syst_upto, size=nsysts)
    nhistograms = [np.random.randint(1, nbins_upto + 1, size=n) for n in nsyst_histos]
    random_alphas = [np.random.uniform(-1, 1, size=nalphas) for n in nsyst_histos]

    random_histogramssets = [
        [  # all histos affected by systematic $nh
            [  # sample $i, systematic $nh
                np.random.uniform(10 * i + j, 10 * i + j + 1, size=nbin).tolist()
                for j in range(3)
            ]
            for i, nbin in enumerate(nh)
        ]
        for nh in nhistograms
    ]
    h, a = filled_shapes(random_histogramssets, random_alphas)
    return h, a

The (slow) interpolations

In all cases, the way we do interpolations is as follows:

1. Loop over both the histogramssets and alphasets simultaneously (e.g. using python’s zip())

2. Loop over all histograms set in the set of histograms sets that correspond to the histograms affected by a given systematic

3. Loop over all of the alphas in the set of alphas

4. Loop over all the bins in the histogram sets simultaneously (e.g. using python’s zip())

5. Apply the interpolation across the same bin index

This is already exhausting to think about, so let’s put this in code form. Depending on the kind of interpolation being done, we’ll pass in func as an argument to the top-level interpolation loop to switch between linear (interpcode=0) and non-linear (interpcode=1).

def interpolation_looper(histogramssets, alphasets, func):
    all_results = []
    for histoset, alphaset in zip(histogramssets, alphasets):
        all_results.append([])
        set_result = all_results[-1]
        for histo in histoset:
            set_result.append([])
            histo_result = set_result[-1]
            for alpha in alphaset:
                alpha_result = []
                for down, nom, up in zip(histo[0], histo[1], histo[2]):
                    v = func(down, nom, up, alpha)
                    alpha_result.append(v)
                histo_result.append(alpha_result)
    return all_results

And we can also define our linear and non-linear interpolations we’ll consider in this notebook that we wish to tensorize.

def interpolation_linear(histogramssets, alphasets):
    def summand(down, nom, up, alpha):
        delta_up = up - nom
        delta_down = nom - down
        if alpha > 0:
            delta = delta_up * alpha
        else:
            delta = delta_down * alpha
        return nom + delta

    return interpolation_looper(histogramssets, alphasets, summand)


def interpolation_nonlinear(histogramssets, alphasets):
    def product(down, nom, up, alpha):
        delta_up = up / nom
        delta_down = down / nom
        if alpha > 0:
            delta = delta_up ** alpha
        else:
            delta = delta_down ** (-alpha)
        return nom * delta

    return interpolation_looper(histogramssets, alphasets, product)

We will also define a helper function that allows us to pass in two functions we wish to compare the outputs for:

def compare_fns(func1, func2):
    h, a = random_histosets_alphasets_pair()

    def _func_runner(func, histssets, alphasets):
        return np.asarray(func(histssets, alphasets))

    old = _func_runner(func1, h, a)
    new = _func_runner(func2, h, a)

    return (np.all(old[~np.isnan(old)] == new[~np.isnan(new)]), (h, a))

For the rest of the notebook, we will detail in explicit form how the linear interpolator gets tensorized, step-by-step. The same sequence of steps will be shown for the non-linear interpolator – but it is left up to the reader to understand the steps.

Tensorizing the Linear Interpolator

Step 0

Step 0 requires converting the innermost conditional check on alpha > 0 into something tensorizable. This also means the calculation itself is going to become tensorized. So we will convert from

if alpha > 0:
    delta =  delta_up*alpha
else:
    delta =  delta_down*alpha

to

delta = np.where(alpha > 0, delta_up*alpha, delta_down*alpha)

Let’s make that change now, and let’s check to make sure we still do the calculation correctly.

# get the internal calculation to use tensorlib backend
def new_interpolation_linear_step0(histogramssets, alphasets):
    all_results = []
    for histoset, alphaset in zip(histogramssets, alphasets):
        all_results.append([])
        set_result = all_results[-1]
        for histo in histoset:
            set_result.append([])
            histo_result = set_result[-1]
            for alpha in alphaset:
                alpha_result = []
                for down, nom, up in zip(histo[0], histo[1], histo[2]):
                    delta_up = up - nom
                    delta_down = nom - down
                    delta = np.where(alpha > 0, delta_up * alpha, delta_down * alpha)
                    v = nom + delta
                    alpha_result.append(v)
                histo_result.append(alpha_result)
    return all_results

And does the calculation still match?

result, (h, a) = compare_fns(interpolation_linear, new_interpolation_linear_step0)
print(result)

True

%%timeit
interpolation_linear(h, a)

189 ms ± 6.14 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%%timeit
new_interpolation_linear_step0(h, a)

255 ms ± 11.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Great! We’re a little bit slower right now, but that’s expected. We’re just getting started.

Step 1

In this step, we would like to remove the innermost zip() call over the histogram bins by calculating the interpolation between the histograms in one fell swoop. This means, instead of writing something like

for down,nom,up in zip(histo[0],histo[1],histo[2]):
    delta_up = up - nom
    ...

one can instead write

delta_up = histo[2] - histo[1]
...

taking advantage of the automatic broadcasting of operations on input tensors. This sort of feature of the tensor backends allows us to speed up code, such as interpolation.

# update the delta variations to remove the zip() call and remove most-nested loop
def new_interpolation_linear_step1(histogramssets, alphasets):
    all_results = []
    for histoset, alphaset in zip(histogramssets, alphasets):
        all_results.append([])
        set_result = all_results[-1]
        for histo in histoset:
            set_result.append([])
            histo_result = set_result[-1]
            for alpha in alphaset:
                alpha_result = []
                deltas_up = histo[2] - histo[1]
                deltas_dn = histo[1] - histo[0]
                calc_deltas = np.where(alpha > 0, deltas_up * alpha, deltas_dn * alpha)
                v = histo[1] + calc_deltas
                alpha_result.append(v)
                histo_result.append(alpha_result)
    return all_results

And does the calculation still match?

result, (h, a) = compare_fns(interpolation_linear, new_interpolation_linear_step1)
print(result)

True

%%timeit
interpolation_linear(h, a)

188 ms ± 7.14 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%%timeit
new_interpolation_linear_step1(h, a)

492 ms ± 42.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Great!

Step 2

In this step, we would like to move the giant array of the deltas calculated to the beginning – outside of all loops – and then only take a subset of it for the calculation itself. This allows us to figure out the entire structure of the input for the rest of the calculations as we slowly move towards including einsum() calls (einstein summation). This means we would like to go from

for histo in histoset:
    delta_up = histo[2] - histo[1]
...

to

all_deltas = ...
for nh, histo in enumerate(histoset):
    deltas = all_deltas[nh]
    ...

Again, we are taking advantage of the automatic broadcasting of operations on input tensors to calculate all the deltas in a single action.

# figure out the giant array of all deltas at the beginning and only take subsets of it for the calculation
def new_interpolation_linear_step2(histogramssets, alphasets):
    all_results = []

    allset_all_histo_deltas_up = histogramssets[:, :, 2] - histogramssets[:, :, 1]
    allset_all_histo_deltas_dn = histogramssets[:, :, 1] - histogramssets[:, :, 0]

    for nset, (histoset, alphaset) in enumerate(zip(histogramssets, alphasets)):
        set_result = []

        all_histo_deltas_up = allset_all_histo_deltas_up[nset]
        all_histo_deltas_dn = allset_all_histo_deltas_dn[nset]

        for nh, histo in enumerate(histoset):
            alpha_deltas = []
            for alpha in alphaset:
                alpha_result = []
                deltas_up = all_histo_deltas_up[nh]
                deltas_dn = all_histo_deltas_dn[nh]
                calc_deltas = np.where(alpha > 0, deltas_up * alpha, deltas_dn * alpha)
                alpha_deltas.append(calc_deltas)
            set_result.append([histo[1] + d for d in alpha_deltas])
        all_results.append(set_result)
    return all_results

And does the calculation still match?

result, (h, a) = compare_fns(interpolation_linear, new_interpolation_linear_step2)
print(result)

True

%%timeit
interpolation_linear(h, a)

179 ms ± 12.4 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_linear_step2(h, a)

409 ms ± 20.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Great!

Step 3

In this step, we get to introduce einstein summation to generalize the calculations we perform across many dimensions in a more concise, straightforward way. See this blog post for some more details on einstein summation notation. In short, it allows us to write

\[c_j = \sum_i \sum_k = A_{ik} B_{kj} \qquad \rightarrow \qquad \texttt{einsum("ij,jk->i", A, B)}\]

in a much more elegant way to express many kinds of common tensor operations such as dot products, transposes, outer products, and so on. This step is generally the hardest as one needs to figure out the corresponding einsum that keeps the calculation preserved (and matching). To some extent it requires a lot of trial and error until you get a feel for how einstein summation notation works.

As a concrete example of a conversion, we wish to go from something like

for nh,histo in enumerate(histoset):
    for alpha in alphaset:
        deltas_up    = all_histo_deltas_up[nh]
        deltas_dn    = all_histo_deltas_dn[nh]
        calc_deltas  = np.where(alpha > 0, deltas_up*alpha, deltas_dn*alpha)
        ...

to get rid of the loop over alpha

for nh,histo in enumerate(histoset):
    alphas_times_deltas_up = np.einsum('i,j->ij',alphaset,all_histo_deltas_up[nh])
    alphas_times_deltas_dn = np.einsum('i,j->ij',alphaset,all_histo_deltas_dn[nh])
    masks = np.einsum('i,j->ij',alphaset > 0,np.ones_like(all_histo_deltas_dn[nh]))

    alpha_deltas  = np.where(masks,alphas_times_deltas_up, alphas_times_deltas_dn)
    ...

In this particular case, we need an outer product that multiplies across the alphaset to the corresponding histoset for the up/down variations. Then we just need to select from either the up variation calculation or the down variation calculation based on the sign of alpha. Try to convince yourself that the einstein summation does what the for-loop does, but a little bit more concisely, and perhaps more clearly! How does the function look now?

# remove the loop over alphas, starts using einsum to help generalize to more dimensions
def new_interpolation_linear_step3(histogramssets, alphasets):
    all_results = []

    allset_all_histo_deltas_up = histogramssets[:, :, 2] - histogramssets[:, :, 1]
    allset_all_histo_deltas_dn = histogramssets[:, :, 1] - histogramssets[:, :, 0]

    for nset, (histoset, alphaset) in enumerate(zip(histogramssets, alphasets)):
        set_result = []

        all_histo_deltas_up = allset_all_histo_deltas_up[nset]
        all_histo_deltas_dn = allset_all_histo_deltas_dn[nset]

        for nh, histo in enumerate(histoset):
            alphas_times_deltas_up = np.einsum(
                'i,j->ij', alphaset, all_histo_deltas_up[nh]
            )
            alphas_times_deltas_dn = np.einsum(
                'i,j->ij', alphaset, all_histo_deltas_dn[nh]
            )
            masks = np.einsum(
                'i,j->ij', alphaset > 0, np.ones_like(all_histo_deltas_dn[nh])
            )

            alpha_deltas = np.where(
                masks, alphas_times_deltas_up, alphas_times_deltas_dn
            )
            set_result.append([histo[1] + d for d in alpha_deltas])

        all_results.append(set_result)
    return all_results

And does the calculation still match?

result, (h, a) = compare_fns(interpolation_linear, new_interpolation_linear_step3)
print(result)

True

%%timeit
interpolation_linear(h, a)

166 ms ± 11.6 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_linear_step3(h, a)

921 ms ± 133 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Great! Note that we’ve been getting a little bit slower during these steps. It will all pay off in the end when we’re fully tensorized! A lot of the internal steps are overkill with the heavy einstein summation and broadcasting at the moment, especially for how many loops in we are.

Step 4

Now in this step, we will move the einstein summations to the outer loop, so that we’re calculating it once! This is the big step, but a little bit easier because all we’re doing is adding extra dimensions into the calculation. The underlying calculation won’t have changed. At this point, we’ll also rename from i and j to a and b for alpha and bin (as in the bin in the histogram). To continue the notation as well, here’s a summary of the dimensions involved:

• s will be for the set under consideration (e.g. the modifier)

• a will be for the alpha variation

• h will be for the histogram affected by the modifier

• b will be for the bin of the histogram

So we wish to move the einsum code from

for nset,(histoset, alphaset) in enumerate(zip(histogramssets,alphasets)):
    ...

    for nh,histo in enumerate(histoset):
        alphas_times_deltas_up = np.einsum('i,j->ij',alphaset,all_histo_deltas_up[nh])
            ...

to

all_alphas_times_deltas_up = np.einsum('...',alphaset,all_histo_deltas_up)
for nset,(histoset, alphaset) in enumerate(zip(histogramssets,alphasets)):
    ...

    for nh,histo in enumerate(histoset):
        ...

So how does this new function look?

# move the einsums to outer loops to get ready to get rid of all loops
def new_interpolation_linear_step4(histogramssets, alphasets):
    allset_all_histo_deltas_up = histogramssets[:, :, 2] - histogramssets[:, :, 1]
    allset_all_histo_deltas_dn = histogramssets[:, :, 1] - histogramssets[:, :, 0]
    allset_all_histo_nom = histogramssets[:, :, 1]

    allsets_all_histos_alphas_times_deltas_up = np.einsum(
        'sa,shb->shab', alphasets, allset_all_histo_deltas_up
    )
    allsets_all_histos_alphas_times_deltas_dn = np.einsum(
        'sa,shb->shab', alphasets, allset_all_histo_deltas_dn
    )
    allsets_all_histos_masks = np.einsum(
        'sa,s...u->s...au', alphasets > 0, np.ones_like(allset_all_histo_deltas_dn)
    )

    allsets_all_histos_deltas = np.where(
        allsets_all_histos_masks,
        allsets_all_histos_alphas_times_deltas_up,
        allsets_all_histos_alphas_times_deltas_dn,
    )

    all_results = []
    for nset, histoset in enumerate(histogramssets):
        all_histos_deltas = allsets_all_histos_deltas[nset]
        set_result = []
        for nh, histo in enumerate(histoset):
            set_result.append([d + histoset[nh, 1] for d in all_histos_deltas[nh]])
        all_results.append(set_result)
    return all_results

And does the calculation still match?

result, (h, a) = compare_fns(interpolation_linear, new_interpolation_linear_step4)
print(result)

True

%%timeit
interpolation_linear(h, a)

160 ms ± 5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_linear_step4(h, a)

119 ms ± 3.19 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Great! And look at that huge speed up in time already, just from moving the multiple, heavy einstein summation calculations up through the loops. We still have some more optimizing to do as we still have explicit loops in our code. Let’s keep at it, we’re almost there!

Step 5

The hard part is mostly over. We have to now think about the nominal variations. Recall that we were trying to add the nominals to the deltas in order to compute the new value. In practice, we’ll return the delta variation only, but we’ll show you how to get rid of this last loop. In this case, we want to figure out how to change code like

all_results    = []
for nset,histoset in enumerate(histogramssets):
    all_histos_deltas = allsets_all_histos_deltas[nset]
    set_result = []
    for nh,histo in enumerate(histoset):
        set_result.append([d + histoset[nh,1] for d in all_histos_deltas[nh]])
    all_results.append(set_result)

to get rid of that most-nested loop

all_results    = []
for nset,histoset in enumerate(histogramssets):
    # look ma, no more loops inside!

So how does this look?

# slowly getting rid of our loops to build the right output tensor -- gotta think about nominals
def new_interpolation_linear_step5(histogramssets, alphasets):
    allset_all_histo_deltas_up = histogramssets[:, :, 2] - histogramssets[:, :, 1]
    allset_all_histo_deltas_dn = histogramssets[:, :, 1] - histogramssets[:, :, 0]
    allset_all_histo_nom = histogramssets[:, :, 1]

    allsets_all_histos_alphas_times_deltas_up = np.einsum(
        'sa,shb->shab', alphasets, allset_all_histo_deltas_up
    )
    allsets_all_histos_alphas_times_deltas_dn = np.einsum(
        'sa,shb->shab', alphasets, allset_all_histo_deltas_dn
    )
    allsets_all_histos_masks = np.einsum(
        'sa,s...u->s...au', alphasets > 0, np.ones_like(allset_all_histo_deltas_dn)
    )

    allsets_all_histos_deltas = np.where(
        allsets_all_histos_masks,
        allsets_all_histos_alphas_times_deltas_up,
        allsets_all_histos_alphas_times_deltas_dn,
    )

    all_results = []

    for nset, (_, alphaset) in enumerate(zip(histogramssets, alphasets)):
        all_histos_deltas = allsets_all_histos_deltas[nset]
        noms = histogramssets[nset, :, 1]

        all_histos_noms_repeated = np.einsum('a,hn->han', np.ones_like(alphaset), noms)

        set_result = all_histos_deltas + all_histos_noms_repeated
        all_results.append(set_result)
    return all_results

And does the calculation still match?

result, (h, a) = compare_fns(interpolation_linear, new_interpolation_linear_step5)
print(result)

True

%%timeit
interpolation_linear(h, a)

160 ms ± 8.28 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_linear_step5(h, a)

1.57 ms ± 75.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Fantastic! And look at the speed up. We’re already faster than the for-loop and we’re not even done yet.

Step 6

The final frontier. Also probably the best Star Wars episode. In any case, we have one more for-loop that needs to die in a slab of carbonite. This should be much easier now that you’re more comfortable with tensor broadcasting and einstein summations.

What does the function look like now?

def new_interpolation_linear_step6(histogramssets, alphasets):
    allset_allhisto_deltas_up = histogramssets[:, :, 2] - histogramssets[:, :, 1]
    allset_allhisto_deltas_dn = histogramssets[:, :, 1] - histogramssets[:, :, 0]
    allset_allhisto_nom = histogramssets[:, :, 1]

    # x is dummy index

    allsets_allhistos_alphas_times_deltas_up = np.einsum(
        'sa,shb->shab', alphasets, allset_allhisto_deltas_up
    )
    allsets_allhistos_alphas_times_deltas_dn = np.einsum(
        'sa,shb->shab', alphasets, allset_allhisto_deltas_dn
    )
    allsets_allhistos_masks = np.einsum(
        'sa,sxu->sxau',
        np.where(alphasets > 0, np.ones(alphasets.shape), np.zeros(alphasets.shape)),
        np.ones(allset_allhisto_deltas_dn.shape),
    )

    allsets_allhistos_deltas = np.where(
        allsets_allhistos_masks,
        allsets_allhistos_alphas_times_deltas_up,
        allsets_allhistos_alphas_times_deltas_dn,
    )
    allsets_allhistos_noms_repeated = np.einsum(
        'sa,shb->shab', np.ones(alphasets.shape), allset_allhisto_nom
    )
    set_results = allsets_allhistos_deltas + allsets_allhistos_noms_repeated
    return set_results

And does the calculation still match?

result, (h, a) = compare_fns(interpolation_linear, new_interpolation_linear_step6)
print(result)

True

%%timeit
interpolation_linear(h, a)

156 ms ± 6.29 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_linear_step6(h, a)

468 µs ± 37.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

And we’re done tensorizing it. There are some more improvements that could be made to make this interpolation calculation even more robust – but for now we’re done.

Tensorizing the Non-Linear Interpolator

This is very, very similar to what we’ve done for the case of the linear interpolator. As such, we will provide the resulting functions for each step, and you can see how things perform all the way at the bottom. Enjoy and learn at your own pace!

def interpolation_nonlinear(histogramssets, alphasets):
    all_results = []
    for histoset, alphaset in zip(histogramssets, alphasets):
        all_results.append([])
        set_result = all_results[-1]
        for histo in histoset:
            set_result.append([])
            histo_result = set_result[-1]
            for alpha in alphaset:
                alpha_result = []
                for down, nom, up in zip(histo[0], histo[1], histo[2]):
                    delta_up = up / nom
                    delta_down = down / nom
                    if alpha > 0:
                        delta = delta_up ** alpha
                    else:
                        delta = delta_down ** (-alpha)
                    v = nom * delta
                    alpha_result.append(v)
                histo_result.append(alpha_result)
    return all_results


def new_interpolation_nonlinear_step0(histogramssets, alphasets):
    all_results = []
    for histoset, alphaset in zip(histogramssets, alphasets):
        all_results.append([])
        set_result = all_results[-1]
        for histo in histoset:
            set_result.append([])
            histo_result = set_result[-1]
            for alpha in alphaset:
                alpha_result = []
                for down, nom, up in zip(histo[0], histo[1], histo[2]):
                    delta_up = up / nom
                    delta_down = down / nom
                    delta = np.where(
                        alpha > 0,
                        np.power(delta_up, alpha),
                        np.power(delta_down, np.abs(alpha)),
                    )
                    v = nom * delta
                    alpha_result.append(v)
                histo_result.append(alpha_result)
    return all_results


def new_interpolation_nonlinear_step1(histogramssets, alphasets):
    all_results = []
    for histoset, alphaset in zip(histogramssets, alphasets):
        all_results.append([])
        set_result = all_results[-1]
        for histo in histoset:
            set_result.append([])
            histo_result = set_result[-1]
            for alpha in alphaset:
                alpha_result = []
                deltas_up = np.divide(histo[2], histo[1])
                deltas_down = np.divide(histo[0], histo[1])
                bases = np.where(alpha > 0, deltas_up, deltas_down)
                exponents = np.abs(alpha)
                calc_deltas = np.power(bases, exponents)
                v = histo[1] * calc_deltas
                alpha_result.append(v)
                histo_result.append(alpha_result)
    return all_results


def new_interpolation_nonlinear_step2(histogramssets, alphasets):
    all_results = []

    allset_all_histo_deltas_up = np.divide(
        histogramssets[:, :, 2], histogramssets[:, :, 1]
    )
    allset_all_histo_deltas_dn = np.divide(
        histogramssets[:, :, 0], histogramssets[:, :, 1]
    )

    for nset, (histoset, alphaset) in enumerate(zip(histogramssets, alphasets)):
        set_result = []

        all_histo_deltas_up = allset_all_histo_deltas_up[nset]
        all_histo_deltas_dn = allset_all_histo_deltas_dn[nset]

        for nh, histo in enumerate(histoset):
            alpha_deltas = []
            for alpha in alphaset:
                alpha_result = []
                deltas_up = all_histo_deltas_up[nh]
                deltas_down = all_histo_deltas_dn[nh]
                bases = np.where(alpha > 0, deltas_up, deltas_down)
                exponents = np.abs(alpha)
                calc_deltas = np.power(bases, exponents)
                alpha_deltas.append(calc_deltas)
            set_result.append([histo[1] * d for d in alpha_deltas])
        all_results.append(set_result)
    return all_results


def new_interpolation_nonlinear_step3(histogramssets, alphasets):
    all_results = []

    allset_all_histo_deltas_up = np.divide(
        histogramssets[:, :, 2], histogramssets[:, :, 1]
    )
    allset_all_histo_deltas_dn = np.divide(
        histogramssets[:, :, 0], histogramssets[:, :, 1]
    )

    for nset, (histoset, alphaset) in enumerate(zip(histogramssets, alphasets)):
        set_result = []

        all_histo_deltas_up = allset_all_histo_deltas_up[nset]
        all_histo_deltas_dn = allset_all_histo_deltas_dn[nset]

        for nh, histo in enumerate(histoset):
            # bases and exponents need to have an outer product, to esentially tile or repeat over rows/cols
            bases_up = np.einsum(
                'a,b->ab', np.ones(alphaset.shape), all_histo_deltas_up[nh]
            )
            bases_dn = np.einsum(
                'a,b->ab', np.ones(alphaset.shape), all_histo_deltas_dn[nh]
            )
            exponents = np.einsum(
                'a,b->ab', np.abs(alphaset), np.ones(all_histo_deltas_up[nh].shape)
            )

            masks = np.einsum(
                'a,b->ab', alphaset > 0, np.ones(all_histo_deltas_dn[nh].shape)
            )
            bases = np.where(masks, bases_up, bases_dn)
            alpha_deltas = np.power(bases, exponents)
            set_result.append([histo[1] * d for d in alpha_deltas])

        all_results.append(set_result)
    return all_results


def new_interpolation_nonlinear_step4(histogramssets, alphasets):
    all_results = []

    allset_all_histo_nom = histogramssets[:, :, 1]
    allset_all_histo_deltas_up = np.divide(
        histogramssets[:, :, 2], allset_all_histo_nom
    )
    allset_all_histo_deltas_dn = np.divide(
        histogramssets[:, :, 0], allset_all_histo_nom
    )

    bases_up = np.einsum(
        'sa,shb->shab', np.ones(alphasets.shape), allset_all_histo_deltas_up
    )
    bases_dn = np.einsum(
        'sa,shb->shab', np.ones(alphasets.shape), allset_all_histo_deltas_dn
    )
    exponents = np.einsum(
        'sa,shb->shab', np.abs(alphasets), np.ones(allset_all_histo_deltas_up.shape)
    )

    masks = np.einsum(
        'sa,shb->shab', alphasets > 0, np.ones(allset_all_histo_deltas_up.shape)
    )
    bases = np.where(masks, bases_up, bases_dn)

    allsets_all_histos_deltas = np.power(bases, exponents)

    all_results = []
    for nset, histoset in enumerate(histogramssets):
        all_histos_deltas = allsets_all_histos_deltas[nset]
        set_result = []
        for nh, histo in enumerate(histoset):
            set_result.append([histoset[nh, 1] * d for d in all_histos_deltas[nh]])
        all_results.append(set_result)
    return all_results


def new_interpolation_nonlinear_step5(histogramssets, alphasets):
    all_results = []

    allset_all_histo_nom = histogramssets[:, :, 1]
    allset_all_histo_deltas_up = np.divide(
        histogramssets[:, :, 2], allset_all_histo_nom
    )
    allset_all_histo_deltas_dn = np.divide(
        histogramssets[:, :, 0], allset_all_histo_nom
    )

    bases_up = np.einsum(
        'sa,shb->shab', np.ones(alphasets.shape), allset_all_histo_deltas_up
    )
    bases_dn = np.einsum(
        'sa,shb->shab', np.ones(alphasets.shape), allset_all_histo_deltas_dn
    )
    exponents = np.einsum(
        'sa,shb->shab', np.abs(alphasets), np.ones(allset_all_histo_deltas_up.shape)
    )

    masks = np.einsum(
        'sa,shb->shab', alphasets > 0, np.ones(allset_all_histo_deltas_up.shape)
    )
    bases = np.where(masks, bases_up, bases_dn)

    allsets_all_histos_deltas = np.power(bases, exponents)

    all_results = []
    for nset, (_, alphaset) in enumerate(zip(histogramssets, alphasets)):
        all_histos_deltas = allsets_all_histos_deltas[nset]
        noms = allset_all_histo_nom[nset]
        all_histos_noms_repeated = np.einsum('a,hn->han', np.ones_like(alphaset), noms)
        set_result = all_histos_deltas * all_histos_noms_repeated
        all_results.append(set_result)
    return all_results


def new_interpolation_nonlinear_step6(histogramssets, alphasets):
    all_results = []

    allset_all_histo_nom = histogramssets[:, :, 1]
    allset_all_histo_deltas_up = np.divide(
        histogramssets[:, :, 2], allset_all_histo_nom
    )
    allset_all_histo_deltas_dn = np.divide(
        histogramssets[:, :, 0], allset_all_histo_nom
    )

    bases_up = np.einsum(
        'sa,shb->shab', np.ones(alphasets.shape), allset_all_histo_deltas_up
    )
    bases_dn = np.einsum(
        'sa,shb->shab', np.ones(alphasets.shape), allset_all_histo_deltas_dn
    )
    exponents = np.einsum(
        'sa,shb->shab', np.abs(alphasets), np.ones(allset_all_histo_deltas_up.shape)
    )

    masks = np.einsum(
        'sa,shb->shab', alphasets > 0, np.ones(allset_all_histo_deltas_up.shape)
    )
    bases = np.where(masks, bases_up, bases_dn)

    allsets_all_histos_deltas = np.power(bases, exponents)
    allsets_allhistos_noms_repeated = np.einsum(
        'sa,shb->shab', np.ones(alphasets.shape), allset_all_histo_nom
    )
    set_results = allsets_all_histos_deltas * allsets_allhistos_noms_repeated
    return set_results

result, (h, a) = compare_fns(interpolation_nonlinear, new_interpolation_nonlinear_step0)
print(result)

True

%%timeit
interpolation_nonlinear(h, a)

149 ms ± 9.45 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_nonlinear_step0(h, a)

527 ms ± 29.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

result, (h, a) = compare_fns(interpolation_nonlinear, new_interpolation_nonlinear_step1)
print(result)

True

%%timeit
interpolation_nonlinear(h, a)

150 ms ± 5.21 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_nonlinear_step1(h, a)

456 ms ± 17.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

result, (h, a) = compare_fns(interpolation_nonlinear, new_interpolation_nonlinear_step2)
print(result)

True

%%timeit
interpolation_nonlinear(h, a)

154 ms ± 4.49 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_nonlinear_step2(h, a)

412 ms ± 31 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

result, (h, a) = compare_fns(interpolation_nonlinear, new_interpolation_nonlinear_step3)
print(result)

True

%%timeit
interpolation_nonlinear(h, a)

145 ms ± 5.15 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_nonlinear_step3(h, a)

1.28 s ± 74.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

result, (h, a) = compare_fns(interpolation_nonlinear, new_interpolation_nonlinear_step4)
print(result)

True

%%timeit
interpolation_nonlinear(h, a)

147 ms ± 8.4 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_nonlinear_step4(h, a)

120 ms ± 3.06 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

result, (h, a) = compare_fns(interpolation_nonlinear, new_interpolation_nonlinear_step5)
print(result)

True

%%timeit
interpolation_nonlinear(h, a)

151 ms ± 5.29 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_nonlinear_step5(h, a)

2.65 ms ± 57.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

result, (h, a) = compare_fns(interpolation_nonlinear, new_interpolation_nonlinear_step6)
print(result)

True

%%timeit
interpolation_nonlinear(h, a)

156 ms ± 3.35 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%%timeit
new_interpolation_nonlinear_step6(h, a)

1.49 ms ± 16 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Empirical Test Statistics

In this notebook we will compute test statistics empirically from pseudo-experiment and establish that they behave as assumed in the asymptotic approximation.

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

np.random.seed(0)
plt.rcParams.update({"font.size": 14})

First, we define the statistical model we will study.

• the signal expected event rate is 10 events

• the background expected event rate is 100 events

• a 10% uncertainty is assigned to the background

The expected event rates are chosen to lie comfortably in the asymptotic regime.

model = pyhf.simplemodels.hepdata_like(
    signal_data=[10.0], bkg_data=[100.0], bkg_uncerts=[10.0]
)

The test statistics based on the profile likelihood described in arXiv:1007.1727 cover scanarios for both

• POI allowed to float to negative values (unbounded; \(\mu \in [-10, 10]\))

• POI constrained to non-negative values (bounded; \(\mu \in [0,10]\))

For consistency, test statistics (\(t_\mu, q_\mu\)) associated with bounded POIs are usually denoted with a tilde (\(\tilde{t}_\mu, \tilde{q}_\mu\)).

We set up the bounds for the fit as follows

unbounded_bounds = model.config.suggested_bounds()
unbounded_bounds[model.config.poi_index] = (-10, 10)

bounded_bounds = model.config.suggested_bounds()

Next we draw some synthetic datasets (also referrerd to as “toys” or pseudo-experiments). We will “throw” 300 toys:

true_poi = 1.0
n_toys = 300
toys = model.make_pdf(pyhf.tensorlib.astensor([true_poi, 1.0])).sample((n_toys,))

In the asymptotic treatment the test statistics are described as a function of the data’s best-fit POI value \(\hat\mu\).

So let’s run some fits so we can plots the empirical test statistics against \(\hat\mu\) to observed the emergence of the asymptotic behavior.

pars = np.asarray(
    [pyhf.infer.mle.fit(toy, model, par_bounds=unbounded_bounds) for toy in toys]
)
fixed_params = model.config.suggested_fixed()

We can now calculate all four test statistics described in arXiv:1007.1727

test_poi = 1.0
tmu = np.asarray(
    [
        pyhf.infer.test_statistics.tmu(
            test_poi,
            toy,
            model,
            init_pars=model.config.suggested_init(),
            par_bounds=unbounded_bounds,
            fixed_params=fixed_params,
        )
        for toy in toys
    ]
)

tmu_tilde = np.asarray(
    [
        pyhf.infer.test_statistics.tmu_tilde(
            test_poi,
            toy,
            model,
            init_pars=model.config.suggested_init(),
            par_bounds=bounded_bounds,
            fixed_params=fixed_params,
        )
        for toy in toys
    ]
)

qmu = np.asarray(
    [
        pyhf.infer.test_statistics.qmu(
            test_poi,
            toy,
            model,
            init_pars=model.config.suggested_init(),
            par_bounds=unbounded_bounds,
            fixed_params=fixed_params,
        )
        for toy in toys
    ]
)

qmu_tilde = np.asarray(
    [
        pyhf.infer.test_statistics.qmu_tilde(
            test_poi,
            toy,
            model,
            init_pars=model.config.suggested_init(),
            par_bounds=bounded_bounds,
            fixed_params=fixed_params,
        )
        for toy in toys
    ]
)

Let’s plot all the test statistics we have computed

muhat = pars[:, model.config.poi_index]
muhat_sigma = np.std(muhat)

We can check the asymptotic assumption that \(\hat{\mu}\) is distributed normally around it’s true value \(\mu' = 1\)

fig, ax = plt.subplots()
fig.set_size_inches(7, 5)

ax.set_xlabel(r"$\hat{\mu}$")
ax.set_ylabel("Density")
ax.set_ylim(top=0.5)

ax.hist(muhat, bins=np.linspace(-4, 5, 31), density=True)
ax.axvline(true_poi, label="true poi", color="black", linestyle="dashed")
ax.axvline(np.mean(muhat), label="empirical mean", color="red", linestyle="dashed")
ax.legend();

Here we define the asymptotic profile likelihood test statistics:

\[t_\mu = \frac{(\mu-\hat\mu)^2}{\sigma^2}\]

\[\begin{split}\tilde{t}_\mu = \begin{cases} t_\mu,\;\text{$\hat{\mu}>0$}\\ t_\mu - t_0,\; \text{else} \end{cases}\end{split}\]

\[\begin{split}q_\mu = \begin{cases} t_\mu,\;\text{$\hat{\mu}<\mu$}\\ 0,\; \text{else} \end{cases}\end{split}\]

\[\begin{split}\tilde{q}_\mu = \begin{cases} \tilde{t}_\mu,\;\text{$\hat{\mu}<\mu$}\\ 0,\; \text{else} \end{cases}\end{split}\]

def tmu_asymp(mutest, muhat, sigma):
    return (mutest - muhat) ** 2 / sigma ** 2


def tmu_tilde_asymp(mutest, muhat, sigma):
    a = tmu_asymp(mutest, muhat, sigma)
    b = tmu_asymp(mutest, muhat, sigma) - tmu_asymp(0.0, muhat, sigma)
    return np.where(muhat > 0, a, b)


def qmu_asymp(mutest, muhat, sigma):
    return np.where(
        muhat < mutest, tmu_asymp(mutest, muhat, sigma), np.zeros_like(muhat)
    )


def qmu_tilde_asymp(mutest, muhat, sigma):
    return np.where(
        muhat < mutest, tmu_tilde_asymp(mutest, muhat, sigma), np.zeros_like(muhat)
    )

And now we can compare them to the empirical values:

muhat_asymp = np.linspace(-3, 5)
fig, axarr = plt.subplots(2, 2)
fig.set_size_inches(14, 10)

labels = [r"$t_{\mu}$", "$\\tilde{t}_{\\mu}$", r"$q_{\mu}$", "$\\tilde{q}_{\\mu}$"]
data = [
    (tmu, tmu_asymp),
    (tmu_tilde, tmu_tilde_asymp),
    (qmu, qmu_asymp),
    (qmu_tilde, qmu_tilde_asymp),
]

for ax, (label, d) in zip(axarr.ravel(), zip(labels, data)):
    empirical, asymp_func = d
    ax.scatter(muhat, empirical, alpha=0.2, label=label)
    ax.plot(
        muhat_asymp,
        asymp_func(1.0, muhat_asymp, muhat_sigma),
        label=f"{label} asymptotic",
        c="r",
    )
    ax.set_xlabel(r"$\hat{\mu}$")
    ax.set_ylabel(f"{label}")
    ax.legend(loc="best")

Using Calculators

One low-level functionality of pyhf when it comes to statistical fits is the idea of a calculator to evaluate with asymptotics or toybased hypothesis testing.

This notebook will introduce very quickly what these calculators are meant to do and how they are used internally in the code. We’ll set up a simple model for demonstration and then show how the calculators come into play.

import numpy as np
import pyhf

np.random.seed(0)

model = pyhf.simplemodels.hepdata_like([6], [9], [3])
data = [9] + model.config.auxdata

The high-level API

If the only thing you are interested in is the hypothesis test result you can just run the high-level API to get it:

CLs_obs, CLs_exp = pyhf.infer.hypotest(1.0, data, model, return_expected_set=True)
print(f'CLs_obs = {CLs_obs}')
print(f'CLs_exp = {CLs_exp}')

CLs_obs = 0.1677886052335611
CLs_exp = [array(0.0159689), array(0.05465771), array(0.16778861), array(0.41863467), array(0.74964133)]

The low-level API

Under the hood, the hypothesis test computes test statistics (such as \(q_\mu, \tilde{q}_\mu\)) and uses calculators in order to assess how likely the computed test statistic value is under various hypotheses. The goal is to provide a consistent API that understands how you wish to perform your hypothesis test.

Let’s look at the asymptotics calculator and then do the same thing for the toybased.

Asymptotics

First, let’s create the calculator for asymptotics using the \(\tilde{q}_\mu\) test statistic.

asymp_calc = pyhf.infer.calculators.AsymptoticCalculator(
    data, model, test_stat='qtilde'
)

Now from this, we want to perform the fit and compute the value of the test statistic from which we can get our \(p\)-values:

teststat = asymp_calc.teststatistic(poi_test=1.0)
print(f'qtilde = {teststat}')

qtilde = 0.0

In addition to this, we can ask the calculator for the distributions of the test statistic for the background-only and signal+background hypotheses:

sb_dist, b_dist = asymp_calc.distributions(poi_test=1.0)

From these distributions, we can ask for the \(p\)-value of the test statistic and use this to calculate the \(\mathrm{CL}_\mathrm{s}\) — a “modified” \(p\)-value.

p_sb = sb_dist.pvalue(teststat)
p_b = b_dist.pvalue(teststat)
p_s = p_sb / p_b

print(f'CL_sb = {p_sb}')
print(f'CL_b = {p_b}')
print(f'CL_s = CL_sb / CL_b = {p_s}')

CL_sb = 0.08389430261678055
CL_b = 0.5
CL_s = CL_sb / CL_b = 0.1677886052335611

In a similar procedure, we can do the same thing for the expected \(\mathrm{CL}_\mathrm{s}\) values as well. We need to get the expected value of the test statistic at each \(\pm\sigma\) and then ask for the expected \(p\)-value associated with each value of the test statistic.

teststat_expected = [b_dist.expected_value(i) for i in [2, 1, 0, -1, -2]]
p_expected = [sb_dist.pvalue(t) / b_dist.pvalue(t) for t in teststat_expected]
p_expected

[0.01596890401598493,
 0.05465770873260968,
 0.1677886052335611,
 0.4186346709326618,
 0.7496413276864433]

However, these sorts of steps are somewhat time-consuming and lengthy, and depending on the calculator chosen, may differ a little bit. The calculator API also serves to harmonize the extraction of the observed \(p\)-values:

p_sb, p_b, p_s = asymp_calc.pvalues(teststat, sb_dist, b_dist)

print(f'CL_sb = {p_sb}')
print(f'CL_b = {p_b}')
print(f'CL_s = CL_sb / CL_b = {p_s}')

CL_sb = 0.08389430261678055
CL_b = 0.5
CL_s = CL_sb / CL_b = 0.1677886052335611

and the expected \(p\)-values:

p_exp_sb, p_exp_b, p_exp_s = asymp_calc.expected_pvalues(sb_dist, b_dist)

print(f'exp. CL_sb = {p_exp_sb}')
print(f'exp. CL_b = {p_exp_b}')
print(f'exp. CL_s = CL_sb / CL_b = {p_exp_s}')

exp. CL_sb = [array(0.00036329), array(0.00867173), array(0.0838943), array(0.35221608), array(0.73258689)]
exp. CL_b = [array(0.02275013), array(0.15865525), array(0.5), array(0.84134475), array(0.97724987)]
exp. CL_s = CL_sb / CL_b = [array(0.0159689), array(0.05465771), array(0.16778861), array(0.41863467), array(0.74964133)]

Toy-Based

The calculator API abstracts away a lot of the differences between various strategies, such that it returns what you want, regardless of whether you choose to perform asymptotics or toy-based testing. It hopefully delivers a simple but powerful API for you!

Let’s create a toy-based calculator and “throw” 500 toys.

toy_calc = pyhf.infer.calculators.ToyCalculator(
    data, model, test_stat='qtilde', ntoys=500
)

Like before, we’ll ask for the test statistic. Unlike the asymptotics case, where we compute the Asimov dataset and perform a series of fits, here we are just evaluating the test statistic for the observed data.

teststat = toy_calc.teststatistic(poi_test=1.0)
print(f'qtilde = {teststat}')

qtilde = 1.902590865638981

inits = model.config.suggested_init()
bounds = model.config.suggested_bounds()
fixeds = model.config.suggested_fixed()
pyhf.infer.test_statistics.qmu_tilde(1.0, data, model, inits, bounds, fixeds)

array(1.90259087)

So now the next thing to do is get our distributions. This is where, in the case of toys, we fit each and every single toy that we’ve randomly sampled from our model.

Note, again, that the API for the calculator is the same as in the asymptotics case.

sb_dist, b_dist = toy_calc.distributions(poi_test=1.0)

From these distributions, we can ask for the \(p\)-value of the test statistic and use this to calculate the \(\mathrm{CL}_\mathrm{s}\).

p_sb, p_b, p_s = toy_calc.pvalues(teststat, sb_dist, b_dist)

print(f'CL_sb = {p_sb}')
print(f'CL_b = {p_b}')
print(f'CL_s = CL_sb / CL_b = {p_s}')

CL_sb = 0.084
CL_b = 0.52
CL_s = CL_sb / CL_b = 0.16153846153846155

In a similar procedure, we can do the same thing for the expected \(\mathrm{CL}_\mathrm{s}\) values as well. We need to get the expected value of the test statistic at each \(\pm\sigma\) and then ask for the expected \(p\)-value associated with each value of the test statistic.

p_exp_sb, p_exp_b, p_exp_s = toy_calc.expected_pvalues(sb_dist, b_dist)

print(f'exp. CL_sb = {p_exp_sb}')
print(f'exp. CL_b = {p_exp_b}')
print(f'exp. CL_s = CL_sb / CL_b = {p_exp_s}')

exp. CL_sb = [array(0.), array(0.008), array(0.084), array(0.318), array(1.)]
exp. CL_b = [array(0.02540926), array(0.17), array(0.52), array(0.846), array(1.)]
exp. CL_s = CL_sb / CL_b = [array(0.), array(0.04594333), array(0.16153846), array(0.37939698), array(1.)]

Examples

Try out in Binder!

Notebooks:

ShapeFactor

import logging
import json
import numpy as np
import matplotlib.pyplot as plt

import pyhf
from pyhf.contrib.viz import brazil

logging.basicConfig(level=logging.INFO)

def prep_data(sourcedata):
    spec = {
        'channels': [
            {
                'name': 'signal',
                'samples': [
                    {
                        'name': 'signal',
                        'data': sourcedata['signal']['bindata']['sig'],
                        'modifiers': [
                            {'name': 'mu', 'type': 'normfactor', 'data': None}
                        ],
                    },
                    {
                        'name': 'bkg1',
                        'data': sourcedata['signal']['bindata']['bkg1'],
                        'modifiers': [
                            {
                                'name': 'coupled_shapefactor',
                                'type': 'shapefactor',
                                'data': None,
                            }
                        ],
                    },
                ],
            },
            {
                'name': 'control',
                'samples': [
                    {
                        'name': 'background',
                        'data': sourcedata['control']['bindata']['bkg1'],
                        'modifiers': [
                            {
                                'name': 'coupled_shapefactor',
                                'type': 'shapefactor',
                                'data': None,
                            }
                        ],
                    }
                ],
            },
        ]
    }
    pdf = pyhf.Model(spec)
    data = []
    for channel in pdf.config.channels:
        data += sourcedata[channel]['bindata']['data']
    data = data + pdf.config.auxdata
    return data, pdf

source = {
    "channels": {
        "signal": {
            "binning": [2, -0.5, 1.5],
            "bindata": {
                "data": [220.0, 230.0],
                "bkg1": [100.0, 70.0],
                "sig": [20.0, 20.0],
            },
        },
        "control": {
            "binning": [2, -0.5, 1.5],
            "bindata": {"data": [200.0, 300.0], "bkg1": [100.0, 100.0]},
        },
    }
}

data, pdf = prep_data(source['channels'])
print(f'data: {data}')

init_pars = pdf.config.suggested_init()
print(f'expected data: {pdf.expected_data(init_pars)}')

par_bounds = pdf.config.suggested_bounds()

INFO:pyhf.pdf:Validating spec against schema: model.json
INFO:pyhf.pdf:adding modifier mu (1 new nuisance parameters)
INFO:pyhf.pdf:adding modifier coupled_shapefactor (2 new nuisance parameters)

data: [200.0, 300.0, 220.0, 230.0]
expected data: [100. 100. 120.  90.]

print(f'initialization parameters: {pdf.config.suggested_init()}')

unconpars = pyhf.infer.mle.fit(data, pdf)
print(f'parameters post unconstrained fit: {unconpars}')

initialization parameters: [1.0, 1.0, 1.0]
parameters post unconstrained fit: [1.00004623 1.99998941 3.00000438]

/srv/conda/envs/notebook/lib/python3.7/site-packages/pyhf/tensor/numpy_backend.py:334: RuntimeWarning: divide by zero encountered in log
  return n * np.log(lam) - lam - gammaln(n + 1.0)

obs_limit, exp_limits, (poi_tests, tests) = pyhf.infer.intervals.upperlimit(
    data, pdf, np.linspace(0, 5, 61), level=0.05, return_results=True
)

/srv/conda/envs/notebook/lib/python3.7/site-packages/pyhf/infer/calculators.py:352: RuntimeWarning: invalid value encountered in double_scalars
  teststat = (qmu - qmu_A) / (2 * self.sqrtqmuA_v)

fig, ax = plt.subplots(figsize=(10, 7))
artists = brazil.plot_results(poi_tests, tests, test_size=0.05, ax=ax)
print(f'expected upper limits: {exp_limits}')
print(f'observed upper limit : {obs_limit}')

expected upper limits: [array(0.74138115), array(0.994935), array(1.38451391), array(1.92899382), array(2.59407668)]
observed upper limit : 2.1945969322493744

XML Import/Export

# NB: python -m pip install pyhf[xmlio]
import pyhf

!ls -lavh ../../../validation/xmlimport_input

total 1752
drwxr-xr-x   7 kratsg  staff   238B Oct 16 22:20 .
drwxr-xr-x  21 kratsg  staff   714B Apr  4 14:26 ..
drwxr-xr-x   6 kratsg  staff   204B Feb 27 17:13 config
drwxr-xr-x   7 kratsg  staff   238B Feb 27 23:41 data
-rw-r--r--   1 kratsg  staff   850K Oct 16 22:20 log
drwxr-xr-x  17 kratsg  staff   578B Nov 15 12:24 results
-rw-r--r--   1 kratsg  staff    21K Oct 16 22:20 scan.pdf

Importing

In order to convert HistFactory XML+ROOT to the pyhf JSON spec for likelihoods, you need to point the command-line interface pyhf xml2json at the top-level XML file. Additionally, as the HistFactory XML specification often uses relative paths, you might need to specify the base directory --basedir from which all other files are located, as specified in the top-level XML. The command will be of the format

pyhf xml2json {top-level XML} --basedir {base directory}

This will print the JSON representation of the XML+ROOT specified. If you wish to store this as a JSON file, you simply need to redirect it

pyhf xml2json {top-level XML} --basedir {base directory} > spec.json

!pyhf xml2json --hide-progress ../../../validation/xmlimport_input/config/example.xml --basedir ../../../validation/xmlimport_input | tee xml_importexport.json

{
    "channels": [
        {
            "name": "channel1",
            "samples": [
                {
                    "data": [
                        20.0,
                        10.0
                    ],
                    "modifiers": [
                        {
                            "data": {
                                "hi": 1.05,
                                "lo": 0.95
                            },
                            "name": "syst1",
                            "type": "normsys"
                        },
                        {
                            "data": null,
                            "name": "SigXsecOverSM",
                            "type": "normfactor"
                        }
                    ],
                    "name": "signal"
                },
                {
                    "data": [
                        100.0,
                        0.0
                    ],
                    "modifiers": [
                        {
                            "data": null,
                            "name": "lumi",
                            "type": "lumi"
                        },
                        {
                            "data": [
                                5.000000074505806,
                                0.0
                            ],
                            "name": "staterror_channel1",
                            "type": "staterror"
                        },
                        {
                            "data": {
                                "hi": 1.05,
                                "lo": 0.95
                            },
                            "name": "syst2",
                            "type": "normsys"
                        }
                    ],
                    "name": "background1"
                },
                {
                    "data": [
                        0.0,
                        100.0
                    ],
                    "modifiers": [
                        {
                            "data": null,
                            "name": "lumi",
                            "type": "lumi"
                        },
                        {
                            "data": [
                                0.0,
                                10.0
                            ],
                            "name": "staterror_channel1",
                            "type": "staterror"
                        },
                        {
                            "data": {
                                "hi": 1.05,
                                "lo": 0.95
                            },
                            "name": "syst3",
                            "type": "normsys"
                        }
                    ],
                    "name": "background2"
                }
            ]
        }
    ],
    "data": {
        "channel1": [
            122.0,
            112.0
        ]
    },
    "toplvl": {
        "measurements": [
            {
                "config": {
                    "parameters": [
                        {
                            "auxdata": [
                                1.0
                            ],
                            "bounds": [
                                [
                                    0.5,
                                    1.5
                                ]
                            ],
                            "fixed": true,
                            "inits": [
                                1.0
                            ],
                            "name": "lumi",
                            "sigmas": [
                                0.1
                            ]
                        },
                        {
                            "fixed": true,
                            "name": "alpha_syst1"
                        }
                    ],
                    "poi": "SigXsecOverSM"
                },
                "name": "GaussExample"
            },
            {
                "config": {
                    "parameters": [
                        {
                            "auxdata": [
                                1.0
                            ],
                            "bounds": [
                                [
                                    0.5,
                                    1.5
                                ]
                            ],
                            "fixed": true,
                            "inits": [
                                1.0
                            ],
                            "name": "lumi",
                            "sigmas": [
                                0.1
                            ]
                        },
                        {
                            "fixed": true,
                            "name": "alpha_syst1"
                        }
                    ],
                    "poi": "SigXsecOverSM"
                },
                "name": "GammaExample"
            },
            {
                "config": {
                    "parameters": [
                        {
                            "auxdata": [
                                1.0
                            ],
                            "bounds": [
                                [
                                    0.5,
                                    1.5
                                ]
                            ],
                            "fixed": true,
                            "inits": [
                                1.0
                            ],
                            "name": "lumi",
                            "sigmas": [
                                0.1
                            ]
                        },
                        {
                            "fixed": true,
                            "name": "alpha_syst1"
                        }
                    ],
                    "poi": "SigXsecOverSM"
                },
                "name": "LogNormExample"
            },
            {
                "config": {
                    "parameters": [
                        {
                            "auxdata": [
                                1.0
                            ],
                            "bounds": [
                                [
                                    0.5,
                                    1.5
                                ]
                            ],
                            "fixed": true,
                            "inits": [
                                1.0
                            ],
                            "name": "lumi",
                            "sigmas": [
                                0.1
                            ]
                        },
                        {
                            "fixed": true,
                            "name": "alpha_syst1"
                        }
                    ],
                    "poi": "SigXsecOverSM"
                },
                "name": "ConstExample"
            }
        ],
        "resultprefix": "./results/example"
    }
}

Exporting

In order to convert the pyhf JSON to the HistFactory XML+ROOT spec for likelihoods, you need to point the command-line interface pyhf json2xml at the JSON file you want to convert. As everything is specified in a single file, there is no need to deal with base directories or looking up additional files. This will produce output XML+ROOT in the --output-dir=./ directory (your current working directory), storing XML configs under --specroot=config and the data file under --dataroot=data. The XML configs are prefixed with --resultprefix=FitConfig by default, so that the top-level XML file will be located at {output dir}/{prefix}.xml. The command will be of the format

pyhf json2xml {JSON spec}

Note that the output directory must already exist.

!mkdir -p output
!pyhf json2xml xml_importexport.json --output-dir output
!ls -lavh output/*

/Users/jovyan/pyhf/src/pyhf/writexml.py:120: RuntimeWarning: invalid value encountered in true_divide
  attrs['HistoName'], np.divide(modifierspec['data'], sampledata).tolist()
-rw-r--r--  1 kratsg  staff   822B Apr  9 09:36 output/FitConfig.xml

output/config:
total 8
drwxr-xr-x  3 kratsg  staff   102B Apr  9 09:36 .
drwxr-xr-x  5 kratsg  staff   170B Apr  9 09:36 ..
-rw-r--r--  1 kratsg  staff   1.0K Apr  9 09:36 FitConfig_channel1.xml

output/data:
total 96
drwxr-xr-x  3 kratsg  staff   102B Apr  9 09:36 .
drwxr-xr-x  5 kratsg  staff   170B Apr  9 09:36 ..
-rw-r--r--  1 kratsg  staff    46K Apr  9 09:36 data.root

!rm xml_importexport.json
!rm -rf output/

import pyhf
import pandas
import numpy as np
import altair as alt

Visualization with Altair

Altair is a python API for generating Vega visuazliation specifications. We demonstracte how to use this to build an interactive chart of pyhf results.

Preparing the data

Altair reads the data as a pandas dataframe, so we create one.

model = pyhf.simplemodels.hepdata_like([7], [20], [5])
data = [25] + model.config.auxdata

muscan = np.linspace(0, 5, 31)
results = [
    pyhf.infer.hypotest(mu, data, model, return_expected_set=True) for mu in muscan
]

data = np.concatenate(
    [
        muscan.reshape(-1, 1),
        np.asarray([r[0] for r in results]).reshape(-1, 1),
        np.asarray([r[1] for r in results]).reshape(-1, 5),
    ],
    axis=1,
)
df = pandas.DataFrame(data, columns=["mu", "obs"] + [f"exp_{i}" for i in range(5)])
df.head()

	mu	obs	exp_0	exp_1	exp_2	exp_3	exp_4
0	0.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000
1	0.166667	0.885208	0.670809	0.771258	0.870322	0.949235	0.989385
2	0.333333	0.795986	0.438838	0.581516	0.743696	0.890881	0.975022
3	0.500000	0.726450	0.279981	0.428500	0.623443	0.825621	0.956105
4	0.666667	0.672216	0.174235	0.308524	0.512383	0.754629	0.931866

Defining the Chart

We need to filled areas for the 1,2 sigma bands and two lines for the expected and observed CLs value. For interactivity we add a hovering label of the observed result

band1 = (
    alt.Chart(df)
    .mark_area(opacity=0.5, color="green")
    .encode(x="mu", y="exp_1", y2="exp_3")
)

band2 = (
    alt.Chart(df)
    .mark_area(opacity=0.5, color="yellow")
    .encode(x="mu", y="exp_0", y2="exp_4")
)

line1 = alt.Chart(df).mark_line(color="black").encode(x="mu", y="obs")

line2 = (
    alt.Chart(df).mark_line(color="black", strokeDash=[5, 5]).encode(x="mu", y="exp_2")
)

nearest = alt.selection_single(
    nearest=True, on="mouseover", fields=["mu"], empty="none"
)


point = (
    alt.Chart(df)
    .mark_point(color="black")
    .encode(x="mu", y="obs", opacity=alt.condition(nearest, alt.value(1), alt.value(0)))
    .add_selection(nearest)
)

text = line1.mark_text(align="left", dx=5, dy=-5).encode(
    text=alt.condition(nearest, "obs", alt.value(" "))
)


band2 + band1 + line1 + line2 + point + text

Hello World, pyhf style

Two bin counting experiment with a background uncertainty

import pyhf

Returning the observed and expected \(\mathrm{CL}_{s}\)

model = pyhf.simplemodels.hepdata_like(
    signal_data=[12.0, 11.0], bkg_data=[50.0, 52.0], bkg_uncerts=[3.0, 7.0]
)
data = [51, 48] + model.config.auxdata
test_mu = 1.0
CLs_obs, CLs_exp = pyhf.infer.hypotest(
    test_mu, data, model, test_stat="qtilde", return_expected=True
)
print(f"Observed: {CLs_obs}, Expected: {CLs_exp}")

Observed: 0.052515541856109765, Expected: 0.06445521290832805

Returning the observed \(\mathrm{CL}_{s}\), \(\mathrm{CL}_{s+b}\), and \(\mathrm{CL}_{b}\)

CLs_obs, p_values = pyhf.infer.hypotest(
    test_mu, data, model, test_stat="qtilde", return_tail_probs=True
)
print(f"Observed CL_s: {CLs_obs}, CL_sb: {p_values[0]}, CL_b: {p_values[1]}")

Observed CL_s: 0.052515541856109765, CL_sb: 0.023324961200974572, CL_b: 0.44415349012077077

A reminder that

\[\mathrm{CL}_{s} = \frac{\mathrm{CL}_{s+b}}{\mathrm{CL}_{b}} = \frac{p_{s+b}}{1-p_{b}}\]

assert CLs_obs == p_values[0] / p_values[1]

Returning the expected \(\mathrm{CL}_{s}\) band values

import numpy as np

CLs_obs, CLs_exp_band = pyhf.infer.hypotest(
    test_mu, data, model, test_stat="qtilde", return_expected_set=True
)
print(f"Observed CL_s: {CLs_obs}\n")
for p_value, n_sigma in enumerate(np.arange(-2, 3)):
    print(
        "Expected CL_s{}: {}".format(
            "      " if n_sigma == 0 else f"({n_sigma} σ)",
            CLs_exp_band[p_value],
        )
    )

Observed CL_s: 0.052515541856109765

Expected CL_s(-2 σ): 0.0026064088679947964
Expected CL_s(-1 σ): 0.013820657528619273
Expected CL_s      : 0.06445521290832805
Expected CL_s(1 σ): 0.23526103626937836
Expected CL_s(2 σ): 0.5730418174887743

Multi-bin Poisson

import logging
import json
import math
import numpy as np
import matplotlib.pyplot as plt

import pyhf
from pyhf import Model, optimizer
from pyhf.simplemodels import hepdata_like
from pyhf.contrib.viz import brazil

from scipy.interpolate import griddata
import scrapbook as sb

def plot_histo(ax, binning, data):
    bin_width = (binning[2] - binning[1]) / binning[0]
    bin_leftedges = np.linspace(binning[1], binning[2], binning[0] + 1)[:-1]
    bin_centers = [le + bin_width / 2.0 for le in bin_leftedges]
    ax.bar(bin_centers, data, 1, alpha=0.5)


def plot_data(ax, binning, data):
    errors = [math.sqrt(d) for d in data]
    bin_width = (binning[2] - binning[1]) / binning[0]
    bin_leftedges = np.linspace(binning[1], binning[2], binning[0] + 1)[:-1]
    bin_centers = [le + bin_width / 2.0 for le in bin_leftedges]
    ax.bar(
        bin_centers,
        data,
        0,
        yerr=errors,
        linewidth=0,
        error_kw=dict(ecolor='k', elinewidth=1),
    )
    ax.scatter(bin_centers, data, c='k')

validation_datadir = '../../validation/data'

source = json.load(open(validation_datadir + '/1bin_example1.json'))
model = hepdata_like(
    source['bindata']['sig'], source['bindata']['bkg'], source['bindata']['bkgerr']
)
data = source['bindata']['data'] + model.config.auxdata

init_pars = model.config.suggested_init()
par_bounds = model.config.suggested_bounds()

obs_limit, exp_limits, (poi_tests, tests) = pyhf.infer.intervals.upperlimit(
    data, model, np.linspace(0, 5, 61), level=0.05, return_results=True
)

/srv/conda/envs/notebook/lib/python3.7/site-packages/pyhf/infer/calculators.py:352: RuntimeWarning: invalid value encountered in double_scalars
  teststat = (qmu - qmu_A) / (2 * self.sqrtqmuA_v)

fig, ax = plt.subplots(figsize=(10, 7))
artists = brazil.plot_results(poi_tests, tests, test_size=0.05, ax=ax)
print(f'expected upper limits: {exp_limits}')
print(f'observed upper limit : {obs_limit}')

expected upper limits: [array(1.07644221), array(1.44922838), array(2.01932904), array(2.83213651), array(3.84750318)]
observed upper limit : 2.381026330918668

source = {
    "binning": [2, -0.5, 1.5],
    "bindata": {
        "data": [120.0, 145.0],
        "bkg": [100.0, 150.0],
        "bkgerr": [15.0, 20.0],
        "sig": [30.0, 45.0],
    },
}


my_observed_counts = source['bindata']['data']

model = hepdata_like(
    source['bindata']['sig'], source['bindata']['bkg'], source['bindata']['bkgerr']
)
data = my_observed_counts + model.config.auxdata


binning = source['binning']

nompars = model.config.suggested_init()


bonly_pars = [x for x in nompars]
bonly_pars[model.config.poi_index] = 0.0
nom_bonly = model.expected_data(bonly_pars, include_auxdata=False)

nom_sb = model.expected_data(nompars, include_auxdata=False)

init_pars = model.config.suggested_init()
par_bounds = model.config.suggested_bounds()

print(init_pars)

bestfit_pars = pyhf.infer.mle.fit(data, model, init_pars, par_bounds)
bestfit_cts = model.expected_data(bestfit_pars, include_auxdata=False)

[1.0, 1.0, 1.0]

f, axarr = plt.subplots(1, 3, sharey=True)
f.set_size_inches(12, 4)

plot_histo(axarr[0], binning, nom_bonly)
plot_data(axarr[0], binning, my_observed_counts)
axarr[0].set_xlim(binning[1:])

plot_histo(axarr[1], binning, nom_sb)
plot_data(axarr[1], binning, my_observed_counts)
axarr[1].set_xlim(binning[1:])

plot_histo(axarr[2], binning, bestfit_cts)
plot_data(axarr[2], binning, my_observed_counts)
axarr[2].set_xlim(binning[1:])

plt.ylim(0, 300);

##  DUMMY 2D thing


def signal(m1, m2):
    massscale = 150.0
    minmass = 100.0
    countscale = 2000

    effective_mass = np.sqrt(m1 ** 2 + m2 ** 2)
    return [countscale * np.exp(-(effective_mass - minmass) / massscale), 0]


def CLs(m1, m2):
    signal_counts = signal(m1, m2)
    pdf = hepdata_like(
        signal_counts, source['bindata']['bkg'], source['bindata']['bkgerr']
    )
    try:
        cls_obs, cls_exp_set = pyhf.infer.hypotest(
            1.0, data, pdf, init_pars, par_bounds, return_expected_set=True
        )
        return cls_obs, cls_exp_set, True
    except AssertionError:
        print(f'fit failed for mass points ({m1}, {m2})')
        return None, None, False

nx, ny = 15, 15
grid = grid_x, grid_y = np.mgrid[
    100 : 1000 : complex(0, nx), 100 : 1000 : complex(0, ny)
]
X = grid.T.reshape(nx * ny, 2)
results = [CLs(m1, m2) for m1, m2 in X]

X = np.array([x for x, (_, _, success) in zip(X, results) if success])
yobs = np.array([obs for obs, exp, success in results if success]).flatten()
yexp = [
    np.array([exp[i] for obs, exp, success in results if success]).flatten()
    for i in range(5)
]

int_obs = griddata(X, yobs, (grid_x, grid_y), method='linear')

int_exp = [griddata(X, yexp[i], (grid_x, grid_y), method='linear') for i in range(5)]

plt.contourf(grid_x, grid_y, int_obs, levels=np.linspace(0, 1))
plt.colorbar()

plt.contour(grid_x, grid_y, int_obs, levels=[0.05], colors='w')
for level in int_exp:
    plt.contour(grid_x, grid_y, level, levels=[0.05], colors='w', linestyles='dashed')

plt.scatter(X[:, 0], X[:, 1], c=yobs, vmin=0, vmax=1);

sb.glue("number_2d_successpoints", len(X))

Data type cannot be displayed: application/scrapbook.scrap.json+json

Multibin Coupled HistoSys

%pylab inline

Populating the interactive namespace from numpy and matplotlib

import logging
import json

import pyhf
from pyhf import Model

logging.basicConfig(level=logging.INFO)

def prep_data(sourcedata):
    spec = {
        "channels": [
            {
                "name": "signal",
                "samples": [
                    {
                        "name": "signal",
                        "data": sourcedata["signal"]["bindata"]["sig"],
                        "modifiers": [
                            {"name": "mu", "type": "normfactor", "data": None}
                        ],
                    },
                    {
                        "name": "bkg1",
                        "data": sourcedata["signal"]["bindata"]["bkg1"],
                        "modifiers": [
                            {
                                "name": "coupled_histosys",
                                "type": "histosys",
                                "data": {
                                    "lo_data": sourcedata["signal"]["bindata"][
                                        "bkg1_dn"
                                    ],
                                    "hi_data": sourcedata["signal"]["bindata"][
                                        "bkg1_up"
                                    ],
                                },
                            }
                        ],
                    },
                    {
                        "name": "bkg2",
                        "data": sourcedata["signal"]["bindata"]["bkg2"],
                        "modifiers": [
                            {
                                "name": "coupled_histosys",
                                "type": "histosys",
                                "data": {
                                    "lo_data": sourcedata["signal"]["bindata"][
                                        "bkg2_dn"
                                    ],
                                    "hi_data": sourcedata["signal"]["bindata"][
                                        "bkg2_up"
                                    ],
                                },
                            }
                        ],
                    },
                ],
            },
            {
                "name": "control",
                "samples": [
                    {
                        "name": "background",
                        "data": sourcedata["control"]["bindata"]["bkg1"],
                        "modifiers": [
                            {
                                "name": "coupled_histosys",
                                "type": "histosys",
                                "data": {
                                    "lo_data": sourcedata["control"]["bindata"][
                                        "bkg1_dn"
                                    ],
                                    "hi_data": sourcedata["control"]["bindata"][
                                        "bkg1_up"
                                    ],
                                },
                            }
                        ],
                    }
                ],
            },
        ]
    }
    pdf = Model(spec)
    data = []
    for c in pdf.spec["channels"]:
        data += sourcedata[c["name"]]["bindata"]["data"]
    data = data + pdf.config.auxdata
    return data, pdf

validation_datadir = "../../validation/data"

source = json.load(open(validation_datadir + "/2bin_2channel_coupledhisto.json"))

data, pdf = prep_data(source["channels"])

print(data)

init_pars = pdf.config.suggested_init()
par_bounds = pdf.config.suggested_bounds()

unconpars = pyhf.infer.mle.fit(data, pdf, init_pars, par_bounds)
print(f"parameters post unconstrained fit: {unconpars}")

conpars = pyhf.infer.mle.fixed_poi_fit(0.0, data, pdf, init_pars, par_bounds)
print(f"parameters post constrained fit: {conpars}")

pdf.expected_data(conpars)

[170.0, 220.0, 110.0, 105.0, 0.0]
parameters post unconstrained fit: [1.53170588e-12 2.21657891e+00]
parameters post constrained fit: [0.         2.21655133]

array([116.08275666, 133.24826999, 183.24826999,  98.08967672,
         2.21655133])

def plot_results(test_mus, cls_obs, cls_exp, poi_tests, test_size=0.05):
    plt.plot(poi_tests, cls_obs, c="k")
    for i, c in zip(range(5), ["grey", "grey", "grey", "grey", "grey"]):
        plt.plot(poi_tests, cls_exp[i], c=c)
    plt.plot(poi_tests, [test_size] * len(test_mus), c="r")
    plt.ylim(0, 1)

def invert_interval(test_mus, cls_obs, cls_exp, test_size=0.05):
    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

poi_tests = np.linspace(0, 5, 61)
tests = [
    pyhf.infer.hypotest(
        poi_test, data, pdf, init_pars, par_bounds, return_expected_set=True
    )
    for poi_test in poi_tests
]
cls_obs = np.array([test[0] for test in tests]).flatten()
cls_exp = [np.array([test[1][i] for test in tests]).flatten() for i in range(5)]

print("\n")
plot_results(poi_tests, cls_obs, cls_exp, poi_tests)
invert_interval(poi_tests, cls_obs, cls_exp)

{'exp': [0.3654675198094938,
  0.4882076670368835,
  0.683262284467055,
  0.9650584704888153,
  1.3142329292131938],
 'obs': 0.3932476110375604}

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

np.random.seed(0)
plt.rcParams.update({"font.size": 14})

Running Monte Carlo simulations (toys)

Finding the (expected) significance can involve costly Monte Carlo calculations (“toys”). The asymptotic approximation described in the paper by Cowan, Cranmer, Gross, Vitells: Asymptotic formulae for likelihood-based tests of new physics [arXiv:1007.1727] provides an alternative to these computationally expensive toy calculations.

This notebook demonstrates a reproduction of one of the key plots in the paper using pyhf.

Counting Experiment

Consider a counting experiment where one observes \(n\) events, following a Poisson distribution with expectation value

\[E[n] = \mu s + b\]

with \(s\) expected signal events and \(b\) expected background events, and signal strength parameter \(\mu\). Follow up in the paper to understand more of the math behind this as the notation is being introduced here. What we will show is the distribution of the (alternative) test statistic \(q_1\) (\(\tilde{q}_1\)) calculated under the assumption of the nominal signal model \((\mu=1)\) for data corresponding to the strength parameter of the background-only \((\mu' = 0)\) and signal+background \((\mu' = 1)\) model hypotheses. For the rest of this notebook, we’ll refer to the background-like model \(\mu'=0\) and the signal-like model \(\mu'=1\).

The first thing we will do is set up the pyhf model with \(s=6\) signal events and \(b=9\) background events (adding a Poisson uncertainty on the background).

signal = 6
background = 9
background_uncertainty = 3
model = pyhf.simplemodels.hepdata_like([signal], [background], [background_uncertainty])
print(f"Channels: {model.config.channels}")
print(f"Samples: {model.config.samples}")
print(f"Parameters: {model.config.parameters}")

Channels: ['singlechannel']
Samples: ['background', 'signal']
Parameters: ['mu', 'uncorr_bkguncrt']

This is a single channel with two samples: signal and background. mu here is the signal strength. Next, we need to define the background-like and signal-like p.d.f.s.

# mu' = 0: background-like
pars_bkg = model.config.suggested_init()
pars_bkg[model.config.poi_index] = 0.0
print(f"Background parameters: {list(zip(model.config.parameters, pars_bkg))}")

# mu' = 1: signal-like
pars_sig = model.config.suggested_init()
pars_sig[model.config.poi_index] = 1.0
print(f"Signal parameters: {list(zip(model.config.parameters, pars_sig))}")

# make the pdfs
pdf_bkg = model.make_pdf(pyhf.tensorlib.astensor(pars_bkg))
pdf_sig = model.make_pdf(pyhf.tensorlib.astensor(pars_sig))

Background parameters: [('mu', 0.0), ('uncorr_bkguncrt', 1.0)]
Signal parameters: [('mu', 1.0), ('uncorr_bkguncrt', 1.0)]

Notice that the parameter of interest, \(\mu'\) is set to zero for background-like models and to one for signal-like models.

Running Toys by Hand

Now that we’ve built our pdfs, we can go ahead and randomly (Monte Carlo) sample them. In this case, we want to “run 10,000 pseudo-experiments” (or “throw toys” as particle physicists would say). This means to draw \(n=10000\) samples from the models:

# note: pdf.sample takes in a "shape" N=(10000,) given the number of samples
n_samples = 10000

# mu' = 0
mc_bkg = pdf_bkg.sample((n_samples,))
# mu' = 1
mc_sig = pdf_sig.sample((n_samples,))

print(mc_bkg.shape)
print(mc_sig.shape)

(10000, 2)
(10000, 2)

You’ll notice that the shape for mc_bkg and mc_sig is not the input shape we passed in (10000,) but rather (10000,2)! Why is that? The HistFactory model is a product of many separate pdfs: Poissons representing the main model, and Gaussians representing the auxiliary measurements. In pyhf, this is represented under the hood as a `Simultaneous <https://scikit-hep.org/pyhf/_generated/pyhf.probability.Simultaneous.html>`__ pdf of the main model and the auxiliary model — hence the second dimension.

We can now calculate the test statistic distributions for \(\tilde{q}_1\) given the background-like and signal-like models. This inference step (running the toys) will take some time:

init_pars = model.config.suggested_init()
par_bounds = model.config.suggested_bounds()
fixed_params = model.config.suggested_fixed()

qtilde_bkg = pyhf.tensorlib.astensor(
    [
        pyhf.infer.test_statistics.qmu_tilde(
            1.0, mc, model, init_pars, par_bounds, fixed_params
        )
        for mc in mc_bkg
    ]
)
qtilde_sig = pyhf.tensorlib.astensor(
    [
        pyhf.infer.test_statistics.qmu_tilde(
            1.0, mc, model, init_pars, par_bounds, fixed_params
        )
        for mc in mc_sig
    ]
)

Running Toys using Calculators

However, as you can see, a lot of this is somewhat cumbersome as you need to carry around two pieces of information: one for background-like and one for signal-like. Instead, pyhf provides a statistics calculator API that both simplifies and harmonizes some of this work for you.

This calculator API allows you to:

• compute a test statistic for the observed data

• provide distributions of that test statistic under various hypotheses

These provided distributions additionally have extra functionality to compute a p-value for the observed test statistic.

We will create a toy-based calculator and evaluate the model \((\mu=1)\) for data simulated under background-like hypothesis \((\mu'=0)\) and under the signal-like hypothesis \((\mu'=1)\). This will compute \(\tilde{q}_1\) for both values of \(\mu'\).

toy_calculator_qtilde = pyhf.infer.utils.create_calculator(
    "toybased",
    model.expected_data(pars_sig),
    model,
    ntoys=n_samples,
    test_stat="qtilde",
)
qtilde_sig, qtilde_bkg = toy_calculator_qtilde.distributions(1.0)

To compute \(q_1\), we just need to alleviate the bounds to allow for \(\mu\) (the parameter of interest) to go below zero. Right now, it is set to the default for normfactor which is [0,10] — a very sensible default most of the time. But if the \(\hat{\mu}\) (the maximum likelihood estimator for \(\mu\)) for our model is truly negative, then we should allow the fit to scan negative \(\mu\) values as well.

qmu_bounds = model.config.suggested_bounds()
print(f"Old bounds: {qmu_bounds}")
qmu_bounds[model.config.poi_index] = (-10, 10)
print(f"New bounds: {qmu_bounds}")

Old bounds: [(0, 10), (1e-10, 10.0)]
New bounds: [(-10, 10), (1e-10, 10.0)]

And then run the toys

toy_calculator_qmu = pyhf.infer.utils.create_calculator(
    "toybased",
    model.expected_data(model.config.suggested_init()),
    model,
    par_bounds=qmu_bounds,
    ntoys=n_samples,
    test_stat="q",
)
qmu_sig, qmu_bkg = toy_calculator_qmu.distributions(1.0)

Signal-like:   0%|          | 0/10000 [00:00<?, ?toy/s]/Users/jovyan/pyhf/src/pyhf/tensor/numpy_backend.py:253: RuntimeWarning: invalid value encountered in log
  return n * np.log(lam) - lam - gammaln(n + 1.0)

Now that we’ve ran the toys, we can make the key plots 🙂.

fig, axes = plt.subplots(nrows=1, ncols=2)
for ax in axes:
    ax.set_xticks(np.arange(0, 10))
ax0, ax1 = axes.flatten()

bins = np.linspace(0, 10, 26)

ax0.hist(
    qmu_sig.samples,
    bins=bins,
    histtype="step",
    density=True,
    label="$f(q_1|1)$ signal-like",
    linewidth=2,
)
ax0.hist(
    qmu_bkg.samples,
    bins=bins,
    histtype="step",
    density=True,
    label="$f(q_1|0)$ background-like",
    linewidth=2,
)
ax0.set_xlabel(r"(a) $q_1$", fontsize=18)
ax0.set_ylabel(r"$f\,(q_1|\mu')$", fontsize=18)
ax0.set_title(r"Test statistic $(q_1)$ distributions")
ax0.legend()

ax1.hist(
    qtilde_sig.samples,
    bins=bins,
    histtype="step",
    density=True,
    label=r"$f(\tilde{q}_1|1)$ signal-like",
    linewidth=2,
)
ax1.hist(
    qtilde_bkg.samples,
    bins=bins,
    histtype="step",
    density=True,
    label=r"$f(\tilde{q}_1|0)$ background-like",
    linewidth=2,
)
ax1.set_xlabel(r"(b) $\tilde{q}_1$", fontsize=18)
ax1.set_ylabel(r"$f\,(\tilde{q}_1|\mu')$", fontsize=18)
ax1.set_title(r"Alternative test statistic $(\tilde{q}_1)$ distributions")
ax1.legend()


plt.setp(axes, xlim=(0, 9), ylim=(1e-3, 2), yscale="log")
fig.set_size_inches(14, 6)
fig.tight_layout(pad=2.0)

from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt

import pyhf
import pyhf.readxml
from pyhf.contrib.viz import brazil

import base64
from IPython.core.display import display, HTML
from ipywidgets import interact, fixed

Binned HEP Statistical Analysis in Python

HistFactory

HistFactory is a popular framework to analyze binned event data and commonly used in High Energy Physics. At its core it is a template for building a statistical model from individual binned distribution (‘Histograms’) and variations on them (‘Systematics’) that represent auxiliary measurements (for example an energy scale of the detector which affects the shape of a distribution)

pyhf

pyhf is a work-in-progress standalone implementation of the HistFactory p.d.f. template and an implementation of the test statistics and asymptotic formulae described in the paper by Cowan, Cranmer, Gross, Vitells: Asymptotic formulae for likelihood-based tests of new physics [arXiv:1007.1727].

Models can be defined using JSON specification, but existing models based on the XML + ROOT file scheme are readable as well.

The Demo

The input data for the statistical analysis was built generated using the containerized workflow engine yadage (see demo from KubeCon 2018 [youtube]). Similarly to Binder this utilizes modern container technology for reproducible science. Below you see the execution graph leading up to the model input data at the bottom.

anim = base64.b64encode(open('workflow.gif', 'rb').read()).decode('ascii')
HTML(f'<img src="data:image/gif;base64,{anim}">')

Read in the Model from XML and ROOT

The ROOT files are read using scikit-hep’s uproot module.

spec = pyhf.readxml.parse('meas.xml', Path.cwd())
workspace = pyhf.Workspace(spec)

From the meas.xml spec, we construct a probability density function (p.d.f). As the model includes systematics, it will be a simultaneous joint p.d.f. of the main model (poisson) and constraint model. The latter is defined by the implied “auxiliary measurements”.

pdf = workspace.model(measurement_name='meas')
data = workspace.data(pdf)
# what is the measurement?
workspace.get_measurement(measurement_name='meas')

{'name': 'meas',
 'config': {'poi': 'SigXsecOverSM',
  'parameters': [{'name': 'lumi',
    'auxdata': [1.0],
    'bounds': [[0.5, 1.5]],
    'inits': [1.0],
    'sigmas': [0.1]},
   {'name': 'SigXsecOverSM',
    'bounds': [[0.0, 3.0]],
    'inits': [1.0],
    'fixed': False}]}}

The p.d.f is build from one data-drived “qcd” (or multijet) estimate and two Monte Carlo-based background samples and is parametrized by five parameters: One parameter of interest SigXsecOverSM and four nuisance parameters that affect the shape of the two Monte Carlo background estimates (both weight-only and shape systematics)

print(f'Samples:\n {workspace.samples}')
print(f'Parameters:\n {workspace.parameters}')

Samples:
 ['mc1', 'mc2', 'qcd', 'signal']
Parameters:
 ['SigXsecOverSM', 'lumi', 'mc1_shape_conv', 'mc1_weight_var1', 'mc2_shape_conv', 'mc2_weight_var1']

par_name_dict = {k: v["slice"].start for k, v in pdf.config.par_map.items()}
all_par_settings = {
    n[0]: tuple(m)
    for n, m in zip(
        sorted(reversed(list(par_name_dict.items())), key=lambda x: x[1]),
        pdf.config.suggested_bounds(),
    )
}
default_par_settings = {n[0]: sum(tuple(m)) / 2.0 for n, m in all_par_settings.items()}

def get_mc_counts(pars):
    deltas, factors = pdf._modifications(pars)
    allsum = pyhf.tensorlib.concatenate(
        deltas + [pyhf.tensorlib.astensor(pdf.nominal_rates)]
    )
    nom_plus_delta = pyhf.tensorlib.sum(allsum, axis=0)
    nom_plus_delta = pyhf.tensorlib.reshape(
        nom_plus_delta, (1,) + pyhf.tensorlib.shape(nom_plus_delta)
    )
    allfac = pyhf.tensorlib.concatenate(factors + [nom_plus_delta])
    return pyhf.tensorlib.product(allfac, axis=0)


animate_plot_pieces = None


def init_plot(fig, ax, par_settings):
    global animate_plot_pieces

    nbins = sum(list(pdf.config.channel_nbins.values()))
    x = np.arange(nbins)
    data = np.zeros(nbins)
    items = []
    for i in [3, 2, 1, 0]:
        items.append(ax.bar(x, data, 1, alpha=1.0))
    animate_plot_pieces = (
        items,
        ax.scatter(x, workspace.data(pdf, with_aux=False), c="k", alpha=1.0, zorder=99),
    )


def animate(ax=None, fig=None, **par_settings):
    global animate_plot_pieces
    items, obs = animate_plot_pieces
    pars = pyhf.tensorlib.astensor(pdf.config.suggested_init())
    for k, v in par_settings.items():
        pars[par_name_dict[k]] = v

    mc_counts = get_mc_counts(pars)
    rectangle_collection = zip(*map(lambda x: x.patches, items))

    for rectangles, binvalues in zip(rectangle_collection, mc_counts[:, 0].T):
        offset = 0
        for sample_index in [3, 2, 1, 0]:
            rect = rectangles[sample_index]
            binvalue = binvalues[sample_index]
            rect.set_y(offset)
            rect.set_height(binvalue)
            offset += rect.get_height()

    fig.canvas.draw()


def plot(ax=None, order=[3, 2, 1, 0], **par_settings):
    pars = pyhf.tensorlib.astensor(pdf.config.suggested_init())
    for k, v in par_settings.items():
        pars[par_name_dict[k]] = v

    mc_counts = get_mc_counts(pars)
    bottom = None
    # nb: bar_data[0] because evaluating only one parset
    for i, sample_index in enumerate(order):
        data = mc_counts[sample_index][0]
        x = np.arange(len(data))
        ax.bar(x, data, 1, bottom=bottom, alpha=1.0)
        bottom = data if i == 0 else bottom + data
    ax.scatter(x, workspace.data(pdf, with_aux=False), c="k", alpha=1.0, zorder=99)

Interactive Exploration of a HistFactory Model

One advantage of a pure-python implementation of Histfactory is the ability to explore the pdf interactively within the setting of a notebook. Try moving the sliders and oberserve the effect on the samples. For example changing the parameter of interest SigXsecOverSM (or µ) controls the overall normalization of the (BSM) signal sample (µ=0 for background-only and µ=1 for the nominal signal-plus-background hypothesis)

%matplotlib notebook
fig, ax = plt.subplots(1, 1)
fig.set_size_inches(10, 5)
ax.set_ylim(0, 1.5 * np.max(workspace.data(pdf, with_aux=False)))

init_plot(fig, ax, default_par_settings)
interact(animate, fig=fixed(fig), ax=fixed(ax), **all_par_settings);

nominal = pdf.config.suggested_init()
background_only = pdf.config.suggested_init()
background_only[pdf.config.poi_index] = 0.0
best_fit = pyhf.infer.mle.fit(data, pdf)

/srv/conda/envs/notebook/lib/python3.7/site-packages/pyhf/tensor/numpy_backend.py:334: RuntimeWarning: invalid value encountered in log
  return n * np.log(lam) - lam - gammaln(n + 1.0)

Fitting

We can now fit the statistical model to the observed data. The best fit of the signal strength is close to the background-only hypothesis.

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, sharex=True)
fig.set_size_inches(18, 4)
ax1.set_ylim(0, 1.5 * np.max(workspace.data(pdf, with_aux=False)))
ax1.set_title('nominal signal + background µ = 1')
plot(ax=ax1, **{k: nominal[v] for k, v in par_name_dict.items()})

ax2.set_title('nominal background-only µ = 0')
plot(ax=ax2, **{k: background_only[v] for k, v in par_name_dict.items()})

ax3.set_title('best fit µ = {:.3g}'.format(best_fit[pdf.config.poi_index]))
plot(ax=ax3, **{k: best_fit[v] for k, v in par_name_dict.items()});

Interval Estimation (Computing Upper Limits on µ)

A common task in the statistical evaluation of High Energy Physics data analyses is the estimation of confidence intervals of parameters of interest. The general strategy is to perform a series of hypothesis tests and then invert the tests in order to obtain an interval with the correct coverage properties.

A common figure of merit is a modified p-value, CLs. Here we compute an upper limit based on a series of CLs tests.

mu_tests = np.linspace(0, 1, 16)
obs_limit, exp_limits, (poi_tests, tests) = pyhf.infer.intervals.upperlimit(
    data, pdf, mu_tests, level=0.05, return_results=True
)

fig, ax = plt.subplots()
fig.set_size_inches(7, 5)

ax.set_title("Hypothesis Tests")
artists = brazil.plot_results(mu_tests, tests, test_size=0.05, ax=ax)

print(f"Observed upper limit: {obs_limit:.3f}\n")
for i, n_sigma in enumerate(np.arange(-2, 3)):
    print(
        "Expected Limit{}: {:.3f}".format(
            "" if n_sigma == 0 else f"({n_sigma} σ)", exp_limits[i]
        )
    )

Observed upper limit: 0.630

Expected Limit(-2 σ): 0.297
Expected Limit(-1 σ): 0.393
Expected Limit: 0.546
Expected Limit(1 σ): 0.762
Expected Limit(2 σ): 1.000

Outreach

We are always interested in talking about pyhf. See the abstract and a list of previously given presentations and feel free to invite us to your next conference/workshop/meeting!

Abstract

The HistFactory p.d.f. template [CERN-OPEN-2012-016] is per-se independent of its implementation in ROOT and it is useful to be able to run statistical analysis outside of the ROOT, RooFit, RooStats framework. pyhf is a pure-python implementation of that statistical model for multi-bin histogram-based analysis and its interval estimation is based on the asymptotic formulas of “Asymptotic formulae for likelihood-based tests of new physics” [1007.1727]. pyhf supports modern computational graph libraries such as TensorFlow and PyTorch in order to make use of features such as auto-differentiation and GPU acceleration.

The HistFactory p.d.f. template
\href{https://cds.cern.ch/record/1456844}{[CERN-OPEN-2012-016]} is
per-se independent of its implementation in ROOT and it is useful to be
able to run statistical analysis outside of the ROOT, RooFit, RooStats
framework. pyhf is a pure-python implementation of that statistical
model for multi-bin histogram-based analysis and its interval
estimation is based on the asymptotic formulas of "Asymptotic formulae
for likelihood-based tests of new physics"
\href{https://arxiv.org/abs/1007.1727}{[arXiv:1007.1727]}. pyhf
supports modern computational graph libraries such as TensorFlow and
PyTorch in order to make use of features such as autodifferentiation
and GPU acceleration.

Presentations

This list will be updated with talks given on pyhf:

• Matthew Feickert. Likelihood preservation and statistical reproduction of searches for new physics. CHEP 2019, Nov 2019. URL: https://indico.cern.ch/event/773049/contributions/3476143/.

• Matthew Feickert. pyhf: pure-Python implementation of HistFactory. PyHEP 2019 Workshop, Oct 2019. URL: https://indico.cern.ch/event/833895/contributions/3577824/.

• Matthew Feickert. pyhf: a pure Python statistical fitting library with tensors and autograd. 19th Python in Science Conference (SciPy 2020), July 2020. URL: http://conference.scipy.org/proceedings/scipy2020/slides.html, doi:10.25080/Majora-342d178e-023.

• Matthew Feickert. pyhf: pure-Python implementation of HistFactory with tensors and automatic differentiation. Tools for High Energy Physics and Cosmology 2020 Workshop, Nov 2020. URL: https://indico.cern.ch/event/955391/contributions/4075505/, doi:10.5281/zenodo.4246056.

• Matthew Feickert, Lukas Heinrich, Giordon Stark, and Kyle Cranmer. pyhf: a pure Python implementation of HistFactory with tensors and autograd. DIANA Meeting - pyhf, October 2018. URL: https://indico.cern.ch/event/759480/.

• Matthew Feickert, Lukas Heinrich, Giordon Stark, and Kyle Cranmer. pyhf: pure-Python implementation of HistFactory models with autograd. 2018 US LHC Users Association Meeting, October 2018. URL: https://indico.fnal.gov/event/17566/session/0/contribution/99.

• Matthew Feickert, Lukas Heinrich, Giordon Stark, and Kyle Cranmer. pyhf: pure-Python implementation of HistFactory models with autograd. (Internal) 3rd ATLAS Machine Learning Workshop, October 2018. URL: https://indico.cern.ch/event/735932/contributions/3136879/.

• Matthew Feickert, Lukas Heinrich, Giordon Stark, and Kyle Cranmer. pyhf: pure-Python implementation of HistFactory models with autograd. (Internal) Joint Machine Learning & Statistics Fora Meeting, September 2018. URL: https://indico.cern.ch/event/757657/contributions/3141134/.

• Lukas Heinrich. Gaussian Process Shape Estimation and Systematics. (Internal) Joint Machine Learning & Statistics Fora Meeting, Dec 2018. URL: https://indico.cern.ch/event/777561/contributions/3234669/.

• Lukas Heinrich. HEP in the Cloud Computing and Open Science Era. EP-IT Data science seminar, Oct 2019. URL: https://indico.cern.ch/event/840837/.

• Lukas Heinrich. Likelihoods associated with statistical fits used in searches for new physics on HEPData and use of RECAST. (Internal) ATLAS Weekly Meeting, Nov 2019. URL: https://indico.cern.ch/event/864395/contributions/3642165/.

• Lukas Heinrich. Traditional inference with machine learning tools. 1st Pan-European Advanced School on Statistics in High Energy Physics, Oct 2019. URL: https://indico.desy.de/indico/event/22731/session/4/contribution/19.

• Lukas Heinrich. pyhf: Full Run-2 ATLAS likelihoods. (Internal) Joint Machine Learning & Statistics Fora Meeting, May 2019. URL: https://indico.cern.ch/event/817483/contributions/3412907/.

• Lukas Heinrich, Matthew Feickert, Giordon Stark, and Kyle Cranmer. pyhf: A standalone HistFactory Implementation. (Re)interpreting the results of new physics searches at the LHC Workshop, May 2018. URL: https://indico.cern.ch/event/702612/contributions/2958658/.

• Giordon Stark. Likelihood Preservation and Reproduction. West Coast LHC Jamboree 2019, Oct 2019. URL: https://indico.cern.ch/event/848030/contributions/3616614/.

• Giordon Stark. New techniques for use of public likelihoods for reinterpretation of search results. 27th International Conference on Supersymmetry and Unification of Fundamental Interactions (SUSY2019), May 2019. URL: https://indico.cern.ch/event/746178/contributions/3396797/.

Tutorials

This list will be updated with tutorials and schools given on pyhf:

• Matthew Feickert. pyhf: Accelerating analyses and preserving likelihoods. PyHEP 2020 Workshop (pyhf v0.4.4), Jul 2020. URL: https://indico.cern.ch/event/882824/contributions/3931292/.

• Lukas Heinrich. Introduction to pyhf. (Internal) ATLAS Induction Day + Software Tutorial (pyhf v0.1.2), Oct 2019. URL: https://indico.cern.ch/event/831761/contributions/3484275/.

• Lukas Heinrich. pyhf tutorial. (Internal) ATLAS Induction Day + Software Tutorial (pyhf v0.4.4), Jul 2020. URL: https://indico.cern.ch/event/892952/contributions/3853306/.

• Giordon Stark. ATLAS Exotics + SUSY Workshop 2020 pyhf Tutorial. ATLAS Exotics + SUSY Workshop 2020, September 2020. URL: https://pyhf.github.io/tutorial-ATLAS-SUSY-Exotics-2020/introduction.html.

Posters

This list will be updated with posters presented on pyhf:

• Matthew Feickert, Lukas Heinrich, Giordon Stark, and Kyle Cranmer. pyhf: a pure Python statistical fitting library for High Energy Physics with tensors and autograd. July 2019. 18th Scientific Computing with Python Conference (SciPy 2019). URL: http://conference.scipy.org/proceedings/scipy2019/slides.html, doi:10.25080/Majora-7ddc1dd1-019.

• Matthew Feickert, Lukas Heinrich, Giordon Stark, and Kyle Cranmer. pyhf: pure Python implementation of HistFactory. November 2019. 24th International Conference on computing in High Energy & Nuclear Physics (CHEP 2019). URL: https://indico.cern.ch/event/773049/contributions/3476180/.

• Lukas Heinrich, Matthew Feickert, Giordon Stark, and Kyle Cranmer. pyhf: auto-differentiable binned statistical models. 19th International Workshop on Advanced Computing and Analysis Techniques in Physics Research (ACAT 2019), March 2019. URL: https://indico.cern.ch/event/708041/contributions/3272095/.

In the Media

This list will be updated with media publications featuring pyhf:

• Katarina Anthony. New open release streamlines interactions with theoretical physicists. ATLAS News, December 2019. URL: https://atlas.cern/updates/atlas-news/new-open-likelihoods.

• Katarina Anthony. New open release allows theorists to explore LHC data in a new way. CERN News, January 2020. URL: https://home.cern/news/news/knowledge-sharing/new-open-release-allows-theorists-explore-lhc-data-new-way.

• Stephanie Melchor. ATLAS releases 'full orchestra' of analysis instruments. Symmetry Magazine, January 2021. URL: https://www.symmetrymagazine.org/article/atlas-releases-full-orchestra-of-analysis-instruments.

Installation

To install, we suggest first setting up a virtual environment

# Python3
python3 -m venv pyhf

and activating it

source pyhf/bin/activate

Install latest stable release from PyPI…

… with NumPy backend

python -m pip install pyhf

… with TensorFlow backend

python -m pip install pyhf[tensorflow]

… with PyTorch backend

python -m pip install pyhf[torch]

… with JAX backend

python -m pip install pyhf[jax]

… with all backends

python -m pip install pyhf[backends]

… with xml import/export functionality

python -m pip install pyhf[xmlio]

Install latest development version from GitHub…

… with NumPy backend

python -m pip install --ignore-installed -U "git+https://github.com/scikit-hep/pyhf.git#egg=pyhf"

… with TensorFlow backend

python -m pip install --ignore-installed -U "git+https://github.com/scikit-hep/pyhf.git#egg=pyhf[tensorflow]"

… with PyTorch backend

python -m pip install --ignore-installed -U "git+https://github.com/scikit-hep/pyhf.git#egg=pyhf[torch]"

… with JAX backend

python -m pip install --ignore-installed -U "git+https://github.com/scikit-hep/pyhf.git#egg=pyhf[jax]"

… with all backends

python -m pip install --ignore-installed -U "git+https://github.com/scikit-hep/pyhf.git#egg=pyhf[backends]"

… with xml import/export functionality

python -m pip install --ignore-installed -U "git+https://github.com/scikit-hep/pyhf.git#egg=pyhf[xmlio]"

Updating pyhf

Rerun the installation command. As the upgrade flag, -U, is used then the libraries will be updated.

Developing

To develop, we suggest using virtual environments together with pip or using pipenv. Once the environment is activated, clone the repo from GitHub

git clone https://github.com/scikit-hep/pyhf.git

and install all necessary packages for development

python -m pip install --ignore-installed -U -e .[complete]

Then setup the Git pre-commit hook for Black by running

pre-commit install

as the rev gets updated through time to track changes of different hooks, simply run

pre-commit autoupdate

to have pre-commit install the new version.

Testing

Data Files

A function-scoped fixture called datadir exists for a given test module which will automatically copy files from the associated test modules data directory into a temporary directory for the given test execution. That is, for example, if a test was defined in test_schema.py, then data files located in test_schema/ will be copied to a temporary directory whose path is made available by the datadir fixture. Therefore, one can do:

def test_patchset(datadir):
    data_file = open(datadir.join("test.txt"))
    ...

which will load the copy of text.txt in the temporary directory. This also works for parameterizations as this will effectively sandbox the file modifications made.

TestPyPI

pyhf tests packaging and distributing by publishing each commit to master to TestPyPI. In addition, installation of the latest test release from TestPyPI can be tested by first installing pyhf normally, to ensure all dependencies are installed from PyPI, and then upgrading pyhf to a dev release from TestPyPI

python -m pip install pyhf
python -m pip install --upgrade --extra-index-url https://test.pypi.org/simple/ --pre pyhf

Note

This adds TestPyPI as an additional package index to search when installing. PyPI will still be the default package index pip will attempt to install from for all dependencies, but if a package has a release on TestPyPI that is a more recent release then the package will be installed from TestPyPI instead. Note that dev releases are considered pre-releases, so 0.1.2 is a “newer” release than 0.1.2.dev3.

Publishing

Publishing to PyPI and TestPyPI is automated through the PyPA’s PyPI publish GitHub Action and the pyhf Tag Creator GitHub Actions workflow. A release can be created from any PR created by a core developer by adding a bumpversion tag to it that corresponds to the release type: major, minor, patch. Once the PR is tagged with the label, the GitHub Actions bot will post a comment with information on the actions it will take once the PR is merged. When the PR has been reviewed, approved, and merged, the Tag Creator workflow will automatically create a new release with bump2version and then deploy the release to PyPI.

Context Files and Archive Metadata

The .zenodo.json and codemeta.json files have the version number automatically updated through bump2version, though their additional metadata should be checked periodically by the dev team (probably every release). The codemeta.json file can be generated automatically from a PyPI install of pyhf using codemetapy

codemetapy --no-extras pyhf > codemeta.json

though the author metadata will still need to be checked and revised by hand. The .zenodo.json is currently generated by hand, so it is worth using codemeta.json as a guide to edit it.

FAQ

Frequently Asked Questions about pyhf and its use.

Questions

Where can I ask questions about pyhf use?

If you have a question about the use of pyhf not covered in the documentation, please ask a question on the GitHub Discussions.

If you believe you have found a bug in pyhf, please report it in the GitHub Issues.

How can I get updates on pyhf?

If you’re interested in getting updates from the pyhf dev team and release announcements you can join the pyhf-announcements mailing list.

Is it possible to set the backend from the CLI?

Yes. Use the --backend option for pyhf cls to specify a tensor backend. The default backend is NumPy. For more information see pyhf cls --help.

Does pyhf support Python 2?

No. Like the rest of the Python community, as of January 2020 the latest releases of pyhf no longer support Python 2. The last release of pyhf that was compatible with Python 2.7 is v0.3.4, which can be installed with

python -m pip install pyhf~=0.3

I only have access to Python 2. How can I use pyhf?

It is recommended that pyhf is used as a standalone step in any analysis, and its environment need not be the same as the rest of the analysis. As Python 2 is not supported it is suggested that you setup a Python 3 runtime on whatever machine you’re using. If you’re using a cluster, talk with your system administrators to get their help in doing so. If you are unable to get a Python 3 runtime, versioned Docker images of pyhf are distributed through Docker Hub.

Once you have Python 3 installed, see the Installation page to get started.

I validated my workspace by comparing pyhf and HistFactory, and while the expected CLs matches, the observed CLs is different. Why is this?

Make sure you’re using the right test statistic (\(q\) or \(\tilde{q}\)) in both situations. In HistFactory, the asymptotics calculator, for example, will do something more involved for the observed CLs if you choose a different test statistic.

I ran validation to compare HistFitter and pyhf, but they don’t match exactly. Why not?

pyhf is validated against HistFactory. HistFitter makes some particular implementation choices that pyhf doesn’t reproduce. Instead of trying to compare pyhf and HistFitter you should try to validate them both against HistFactory.

How is pyhf typeset?

As you may have guessed from this page, pyhf is typeset in all lowercase. This is largely historical, as the core developers had just always typed it that way and it seemed a bit too short of a library name to write as PyHF. When typesetting in LaTeX the developers recommend introducing the command

\newcommand{\pyhf}{\texttt{pyhf}}

If the journal you are publishing in requires you to use textsc for software names it is okay to instead use

\newcommand{\pyhf}{\textsc{pyhf}}

Troubleshooting

• import torch or import pyhf causes a Segmentation fault (core dumped)

  This is may be the result of a conflict with the NVIDIA drivers that you have installed on your machine. Try uninstalling and completely removing all of them from your machine

  # On Ubuntu/Debian
  sudo apt-get purge nvidia*
  

  and then installing the latest versions.

Translations

One key goal of pyhf is to provide seamless translations between other statistical frameworks and pyhf. This page details the various ways to translate from a tool you might already be using as part of an existing analysis to pyhf. Many of these solutions involve extracting out the HistFactory workspace and then running pyhf xml2json which provides a single JSON workspace that can be loaded directly into pyhf.

HistFitter

In order to go from HistFitter to pyhf, the first step is to extract out the HistFactory workspaces. Assuming you have an existing configuration file, config.py, you likely run an exclusion fit like so:

HistFitter.py -f -D "before,after,corrMatrix" -F excl config.py

The name of output workspace files depends on four parameters you define in your config.py:

• analysisName is from configMgr.analysisName

• prefix is defined in configMgr.addFitConfig({prefix})

• measurementName is the first measurement you define via fitConfig.addMeasurement(name={measurementName},...)

• channelName are the names of channels you define via fitConfig.addChannel("cuts", [{channelName}], ...)

• cachePath is where HistFitter stores the cached histograms, defined by configMgr.histCacheFile which defaults to data/{analysisName}.root

To dump the HistFactory workspace, you will modify the above to skip the fit -f and plotting -D so you end up with

HistFitter.py -wx -F excl config.py

The -w flag tells HistFitter to (re)create the HistFactory workspace stored in results/{analysisName}/{prefix}_combined_{measurementName}.root. The -x flag tells HistFitter to dump the XML files into config/{analysisName}/, with the top-level file being {prefix}.xml and all other files being {prefix}_{channelName}_cuts.xml.

Typically, prefix = 'FitConfig' and measurementName = 'NormalMeasurement'. For example, if the following exists in your config.py

from configManager import configMgr

# ...
configMgr.analysisName = "3b_tag21.2.27-1_RW_ExpSyst_36100_multibin_bkg"
configMgr.histCacheFile = f"cache/{configMgr.analysisName:s}.root"
# ...
fitConfig = configMgr.addFitConfig("Excl")
# ...
channel = fitConfig.addChannel("cuts", ["SR_0L"], 1, 0.5, 1.5)
# ...
meas1 = fitConfig.addMeasurement(name="DefaultMeasurement", lumi=1.0, lumiErr=0.029)
meas1.addPOI("mu_SIG1")
# ...
meas2 = fitConfig.addMeasurement(name="DefaultMeasurement", lumi=1.0, lumiErr=0.029)
meas2.addPOI("mu_SIG2")

Then, you expect the following files to be made:

• config/3b_tag21.2.27-1_RW_ExpSyst_36100_multibin_bkg/Excl.xml

• config/3b_tag21.2.27-1_RW_ExpSyst_36100_multibin_bkg/Excl_SR_0L_cuts.xml

• cache/3b_tag21.2.27-1_RW_ExpSyst_36100_multibin_bkg.root

• results/3b_tag21.2.27-1_RW_ExpSyst_36100_multibin_bkg/Excl_combined_DefaultMeasurement.root

These are all the files you need in order to use pyhf xml2json. At this point, you could run

pyhf xml2json config/3b_tag21.2.27-1_RW_ExpSyst_36100_multibin_bkg/Excl.xml

which will read all of the XML files and load the histogram data from the histogram cache.

The HistFactory workspace in results/ contains all of the information necessary to rebuild the XML files again. For debugging purposes, the pyhf developers will often ask for your workspace file, which means results/3b_tag21.2.27-1_RW_ExpSyst_36100_multibin_bkg/Excl_combined_DefaultMeasurement.root. If you want to generate the XML, you can open this file in ROOT and run DefaultMeasurement->PrintXML() which puts all of the XML files into the current directory you are in.

TRExFitter

Note

For more details on this section, please refer to the ATLAS-internal TRExFitter documentation.

In order to go from TRExFitter to pyhf, the good news is that the RooWorkspace files (XML and ROOT) are already made for you. For a given configuration which looks like

Job: "pyhf_example"
Label: "..."

You can expect some files to be made after the n/h and w steps:

• pyhf_example/RooStats/pyhf_example.xml

• pyhf_example/RooStats/pyhf_example_Signal_region.xml

• pyhf_example/Histograms/pyhf_example_Signal_region_histos.root

These are all the files you need in order to use pyhf xml2json. At this point, you could run

pyhf xml2json pyhf_example/RooStats/pyhf_example.xml

which will read all of the XML files and load the histogram data from the histogram cache.

Warning

There are a few caveats one needs to be aware of with this conversion:

• Uncorrelated shape systematics cannot be pruned, see Issue #662.

• Custom expressions for normalization factors cannot be used, see Issue #850.

Command Line API

pyhf

Top-level CLI entrypoint.

pyhf [OPTIONS] COMMAND [ARGS]...

Options

--version

Show the version and exit.

--cite, --citation

Print the bibtex citation for this software

cls

Compute CLs value(s) for a given pyhf workspace.

Example:

$ curl -sL https://git.io/JJYDE | pyhf cls


{
    "CLs_exp": [
        0.07807427911686156,
        0.17472571775474618,
        0.35998495263681285,
        0.6343568235898907,
        0.8809947004472013
    ],
    "CLs_obs": 0.3599845631401915
}

pyhf cls [OPTIONS] [WORKSPACE]

Options

--output-file <output_file>

The location of the output json file. If not specified, prints to screen.

--measurement <measurement>

-p, --patch <patch>

--test-poi <test_poi>

--test-stat <test_stat>

Options

q | qtilde

--calctype <calctype>

Options

asymptotics | toybased

--backend <backend>

The tensor backend used for the calculation.

Options

numpy | pytorch | tensorflow | jax | np | torch | tf

--optimizer <optimizer>

The optimizer used for the calculation.

Options

scipy | minuit

--optconf <optconf>

Arguments

WORKSPACE

Optional argument

combine

Combine two workspaces into a single workspace.

See pyhf.workspace.Workspace.combine() for more information.

pyhf combine [OPTIONS] [WORKSPACE_ONE] [WORKSPACE_TWO]

Options

-j, --join <join>

The join operation to apply when combining the two workspaces.

Options

none | outer | left outer | right outer

--output-file <output_file>

The location of the output json file. If not specified, prints to screen.

--merge-channels, --no-merge-channels

Whether or not to deeply merge channels. Can only be done with left/right outer joins.

Arguments

WORKSPACE_ONE

Optional argument

WORKSPACE_TWO

Optional argument

completions

Generate shell completion code.

pyhf completions [OPTIONS] SHELL

Arguments

SHELL

Required argument

contrib

Contrib experimental operations.

Note

Requires installation of the contrib extra.

$ python -m pip install pyhf[contrib]

pyhf contrib [OPTIONS] COMMAND [ARGS]...

download

Download the patchset archive from the remote URL and extract it in a directory at the path given.

Example:

$ pyhf contrib download --verbose https://doi.org/10.17182/hepdata.90607.v3/r3 1Lbb-likelihoods


1Lbb-likelihoods/patchset.json
1Lbb-likelihoods/README.md
1Lbb-likelihoods/BkgOnly.json

Raises:

InvalidArchiveHost: if the provided archive host name is not known to be valid

pyhf contrib download [OPTIONS] ARCHIVE_URL OUTPUT_DIRECTORY

Options

-v, --verbose

Enables verbose mode

-f, --force

Force download from non-approved host

-c, --compress

Keep the archive in a compressed tar.gz form

Arguments

ARCHIVE_URL

Required argument

OUTPUT_DIRECTORY

Required argument

digest

Use hashing algorithm to calculate the workspace digest.

Returns:

digests (dict): A mapping of the hashing algorithms used to the computed digest for the workspace.

Example:

$ curl -sL https://raw.githubusercontent.com/scikit-hep/pyhf/master/docs/examples/json/2-bin_1-channel.json | pyhf digest
sha256:dad8822af55205d60152cbe4303929042dbd9d4839012e055e7c6b6459d68d73

pyhf digest [OPTIONS] [WORKSPACE]

Options

-a, --algorithm <algorithm>

The hashing algorithm used to compute the workspace digest.

-j, --json, -p, --plaintext

Output the hash values as a JSON dictionary or plaintext strings

Arguments

WORKSPACE

Optional argument

fit

Perform a maximum likelihood fit for a given pyhf workspace.

Example:

$ curl -sL https://git.io/JJYDE | pyhf fit --value


{
    "mle_parameters": {
        "mu": [
            0.00017298628839781602
        ],
        "uncorr_bkguncrt": [
            1.0000015671710816,
            0.9999665895859197
        ]
    },
    "twice_nll": 23.19636590468879
}

pyhf fit [OPTIONS] [WORKSPACE]

Options

--output-file <output_file>

The location of the output json file. If not specified, prints to screen.

--measurement <measurement>

-p, --patch <patch>

--value

Flag for returning the fitted value of the objective function.

--backend <backend>

The tensor backend used for the calculation.

Options

numpy | pytorch | tensorflow | jax | np | torch | tf

--optimizer <optimizer>

The optimizer used for the calculation.

Options

scipy | minuit

--optconf <optconf>

Arguments

WORKSPACE

Optional argument

inspect

Inspect a pyhf JSON document.

Example:

$ curl -sL https://raw.githubusercontent.com/scikit-hep/pyhf/master/docs/examples/json/2-bin_1-channel.json | pyhf inspect
          Summary
    ------------------
       channels  1
        samples  2
     parameters  2
      modifiers  2

       channels  nbins
     ----------  -----
  singlechannel    2

        samples
     ----------
     background
         signal

     parameters  constraint              modifiers
     ----------  ----------              ----------
             mu  unconstrained           normfactor
uncorr_bkguncrt  constrained_by_poisson  shapesys

    measurement           poi            parameters
     ----------        ----------        ----------
(*) Measurement            mu            (none)

pyhf inspect [OPTIONS] [WORKSPACE]

Options

--output-file <output_file>

The location of the output json file. If not specified, prints to screen.

--measurement <measurement>

Arguments

WORKSPACE

Optional argument

json2xml

Convert pyhf JSON back to XML + ROOT files.

pyhf json2xml [OPTIONS] [WORKSPACE]

Options

--output-dir <output_dir>

--specroot <specroot>

--dataroot <dataroot>

--resultprefix <resultprefix>

-p, --patch <patch>

Arguments

WORKSPACE

Optional argument

patchset

Operations involving patchsets.

pyhf patchset [OPTIONS] COMMAND [ARGS]...

apply

Apply a patch from patchset to the background-only workspace specification.

Raises:

InvalidPatchLookup: if the provided patch name is not in the patchset PatchSetVerificationError: if the patchset cannot be verified against the workspace specification

Returns:

workspace (Workspace): The patched background-only workspace.

pyhf patchset apply [OPTIONS] [BACKGROUND_ONLY] [PATCHSET]

Options

--name <name>

The name of the patch to extract.

--output-file <output_file>

The location of the output json file. If not specified, prints to screen.

Arguments

BACKGROUND_ONLY

Optional argument

PATCHSET

Optional argument

extract

Extract a patch from a patchset.

Raises:

InvalidPatchLookup: if the provided patch name is not in the patchset

Returns:

jsonpatch (list): A list of jsonpatch operations to apply to a workspace.

pyhf patchset extract [OPTIONS] [PATCHSET]

Options

--name <name>

The name of the patch to extract.

--output-file <output_file>

The location of the output json file. If not specified, prints to screen.

--with-metadata, --without-metadata

Include patchset metadata in output.

Arguments

PATCHSET

Optional argument

inspect

Inspect the PatchSet (e.g. list individual patches).

Returns:

None

pyhf patchset inspect [OPTIONS] [PATCHSET]

Arguments

PATCHSET

Optional argument

verify

Verify the patchset digests against a background-only workspace specification. Verified if no exception was raised.

Raises:

PatchSetVerificationError: if the patchset cannot be verified against the workspace specification

Returns:

None

pyhf patchset verify [OPTIONS] [BACKGROUND_ONLY] [PATCHSET]

Arguments

BACKGROUND_ONLY

Optional argument

PATCHSET

Optional argument

prune

Prune components from the workspace.

See pyhf.workspace.Workspace.prune() for more information.

pyhf prune [OPTIONS] [WORKSPACE]

Options

--output-file <output_file>

The location of the output json file. If not specified, prints to screen.

-c, --channel <CHANNEL>...>

-s, --sample <SAMPLE>...>

-m, --modifier <MODIFIER>...>

-t, --modifier-type <modifier_type>

Options

histosys | lumi | normfactor | normsys | shapefactor | shapesys | staterror

--measurement <MEASUREMENT>...>

Arguments

WORKSPACE

Optional argument

rename

Rename components of the workspace.

See pyhf.workspace.Workspace.rename() for more information.

pyhf rename [OPTIONS] [WORKSPACE]

Options

--output-file <output_file>

The location of the output json file. If not specified, prints to screen.

-c, --channel <PATTERN> <REPLACE>...>

-s, --sample <PATTERN> <REPLACE>...>

-m, --modifier <PATTERN> <REPLACE>...>

--measurement <PATTERN> <REPLACE>...>

Arguments

WORKSPACE

Optional argument

sort

Sort the workspace.

See pyhf.workspace.Workspace.sorted() for more information.

Example:

$ curl -sL https://raw.githubusercontent.com/scikit-hep/pyhf/master/docs/examples/json/2-bin_1-channel.json | pyhf sort | jq '.' | md5
8be5186ec249d2704e14dd29ef05ffb0

$ curl -sL https://raw.githubusercontent.com/scikit-hep/pyhf/master/docs/examples/json/2-bin_1-channel.json | jq -S '.channels|=sort_by(.name)|.channels[].samples|=sort_by(.name)|.channels[].samples[].modifiers|=sort_by(.name,.type)|.observations|=sort_by(.name)' | md5
8be5186ec249d2704e14dd29ef05ffb0

pyhf sort [OPTIONS] [WORKSPACE]

Options

--output-file <output_file>

The location of the output json file. If not specified, prints to screen.

Arguments

WORKSPACE

Optional argument

xml2json

Entrypoint XML: The top-level XML file for the PDF definition.

pyhf xml2json [OPTIONS] ENTRYPOINT_XML

Options

--basedir <basedir>

The base directory for the XML files to point relative to.

--output-file <output_file>

The location of the output json file. If not specified, prints to screen.

--track-progress, --hide-progress

Arguments

ENTRYPOINT_XML

Required argument

Python API

Top-Level

default_backend	NumPy backend for pyhf
default_optimizer	Optimizer that uses scipy.optimize.minimize().
tensorlib	NumPy backend for pyhf
optimizer	Optimizer that uses scipy.optimize.minimize().
get_backend	Get the current backend and the associated optimizer
set_backend	Set the backend and the associated optimizer
readxml
writexml

Probability Distribution Functions (PDFs)

Normal	The Normal distribution with mean loc and standard deviation scale.
Poisson	The Poisson distribution with rate parameter rate.
Independent	A probability density corresponding to the joint distribution of a batch of identically distributed random variables.
Simultaneous	A probability density corresponding to the joint distribution of multiple non-identical component distributions

Making Models from PDFs

Model	The main pyhf model class.
_ModelConfig
Workspace	A JSON-serializable object that is built from an object that follows the workspace.json schema.
PatchSet	A way to store a collection of patches (Patch).
Patch	A way to store a patch definition as part of a patchset (PatchSet).
simplemodels.hepdata_like	Construct a simple single channel Model with a shapesys modifier representing an uncorrelated background uncertainty.

Backends

The computational backends that pyhf provides interfacing for the vector-based calculations.

numpy_backend.numpy_backend	NumPy backend for pyhf
pytorch_backend.pytorch_backend	PyTorch backend for pyhf
tensorflow_backend.tensorflow_backend	TensorFlow backend for pyhf
jax_backend.jax_backend	JAX backend for pyhf

Optimizers

mixins.OptimizerMixin	Mixin Class to build optimizers.
opt_scipy.scipy_optimizer	Optimizer that uses scipy.optimize.minimize().
opt_minuit.minuit_optimizer	Optimizer that minimizes via iminuit.Minuit.migrad().

Modifiers

histosys
normfactor
normsys
shapefactor
shapesys
staterror

Interpolators

code0	The piecewise-linear interpolation strategy.
code1	The piecewise-exponential interpolation strategy.
code2	The quadratic interpolation and linear extrapolation strategy.
code4	The polynomial interpolation and exponential extrapolation strategy.
code4p	The piecewise-linear interpolation strategy, with polynomial at \(\left|a\right| < 1\).

Inference

Test Statistics

test_statistics.q0	The test statistic, \(q_{0}\), for discovery of a positive signal as defined in Equation (12) in [1007.1727], for \(\mu=0\).
test_statistics.qmu	The test statistic, \(q_{\mu}\), for establishing an upper limit on the strength parameter, \(\mu\), as defiend in Equation (14) in [1007.1727]
test_statistics.qmu_tilde	The “alternative” test statistic, \(\tilde{q}_{\mu}\), for establishing an upper limit on the strength parameter, \(\mu\), for models with bounded POI, as defiend in Equation (16) in [1007.1727]
test_statistics.tmu	The test statistic, \(t_{\mu}\), for establishing a two-sided interval on the strength parameter, \(\mu\), as defiend in Equation (8) in [1007.1727]
test_statistics.tmu_tilde	The test statistic, \(\tilde{t}_{\mu}\), for establishing a two-sided interval on the strength parameter, \(\mu\), for models with bounded POI, as defiend in Equation (11) in [1007.1727]
utils.get_test_stat	Get the test statistic function by name.

Calculators

calculators.generate_asimov_data	Compute Asimov Dataset (expected yields at best-fit values) for a given POI value.
calculators.AsymptoticTestStatDistribution	The distribution the test statistic in the asymptotic case.
calculators.EmpiricalDistribution	The empirical distribution of the test statistic.
calculators.AsymptoticCalculator	The Asymptotic Calculator.
calculators.ToyCalculator	The Toy-based Calculator.
utils.create_calculator	Creates a calculator object of the specified calctype.

Fits and Tests

mle.twice_nll	Two times the negative log-likelihood of the model parameters, \(\left(\mu, \boldsymbol{\theta}\right)\), given the observed data.
mle.fit	Run a maximum likelihood fit.
mle.fixed_poi_fit	Run a maximum likelihood fit with the POI value fixed.
hypotest	Compute \(p\)-values and test statistics for a single value of the parameter of interest.
intervals.upperlimit	Calculate an upper limit interval (0, poi_up) for a single Parameter of Interest (POI) using a fixed scan through POI-space.
utils.all_pois_floating	Check whether all POI(s) are floating (i.e.

Exceptions

Various exceptions, apart from standard python exceptions, that are raised from using the pyhf API.

InvalidMeasurement	InvalidMeasurement is raised when a specified measurement is invalid given the specification.
InvalidNameReuse
InvalidSpecification	InvalidSpecification is raised when a specification does not validate against the given schema.
InvalidPatchSet	InvalidPatchSet is raised when a given patchset object does not have the right configuration, even though it validates correctly against the schema.
InvalidPatchLookup	InvalidPatchLookup is raised when the patch lookup from a patchset object has failed
PatchSetVerificationError	PatchSetVerificationError is raised when the workspace digest does not match the patchset digests as part of the verification procedure
InvalidWorkspaceOperation	InvalidWorkspaceOperation is raised when an operation on a workspace fails.
InvalidModel	InvalidModel is raised when a given model does not have the right configuration, even though it validates correctly against the schema.
InvalidModifier	InvalidModifier is raised when an invalid modifier is requested.
InvalidInterpCode	InvalidInterpCode is raised when an invalid/unimplemented interpolation code is requested.
ImportBackendError	MissingLibraries is raised when something is imported by sustained an import error due to missing additional, non-default libraries.
InvalidBackend	InvalidBackend is raised when trying to set a backend that does not exist.
InvalidOptimizer	InvalidOptimizer is raised when trying to set an optimizer that does not exist.
InvalidPdfParameters	InvalidPdfParameters is raised when trying to evaluate a pdf with invalid parameters.
InvalidPdfData	InvalidPdfData is raised when trying to evaluate a pdf with invalid data.

Utilities

load_schema
validate
options_from_eqdelimstring
digest	Get the digest for the provided object.

Contrib

viz.brazil	Brazil Band Plots.
utils.download	Download the patchset archive from the remote URL and extract it in a directory at the path given.

Use and Citations

Warning: This is a development version and should not be cited. To find the specific version to cite, please go to ReadTheDocs.

Citation

The preferred BibTeX entry for citation of pyhf includes both the Zenodo archive and the JOSS paper:

@software{pyhf,
  author = {Lukas Heinrich and Matthew Feickert and Giordon Stark},
  title = "{pyhf: v0.6.1}",
  version = {0.6.1},
  doi = {10.5281/zenodo.1169739},
  url = {https://github.com/scikit-hep/pyhf},
}

@article{pyhf_joss,
  doi = {10.21105/joss.02823},
  url = {https://doi.org/10.21105/joss.02823},
  year = {2021},
  publisher = {The Open Journal},
  volume = {6},
  number = {58},
  pages = {2823},
  author = {Lukas Heinrich and Matthew Feickert and Giordon Stark and Kyle Cranmer},
  title = {pyhf: pure-Python implementation of HistFactory statistical models},
  journal = {Journal of Open Source Software}
}

Use in Publications

Updating list of citations and use cases of pyhf:

• Waleed Abdallah and others. Reinterpretation of LHC Results for New Physics: Status and Recommendations after Run 2. 2020. arXiv:2003.07868.

• Gaël Alguero, Jan Heisig, Charanjit K. Khosa, Sabine Kraml, Suchita Kulkarni, Andre Lessa, Philipp Neuhuber, Humberto Reyes-González, Wolfgang Waltenberger, and Alicia Wongel. New developments in SModelS. In Tools for High Energy Physics and Cosmology. 12 2020. arXiv:2012.08192.

• Gaël Alguero, Sabine Kraml, and Wolfgang Waltenberger. A SModelS interface for pyhf likelihoods. Sep 2020. arXiv:2009.01809.

• J. Alison and others. Higgs Boson Pair Production at Colliders: Status and Perspectives. In B. Di Micco, M. Gouzevitch, J. Mazzitelli, and C. Vernieri, editors, Double Higgs Production at Colliders. 9 2019. arXiv:1910.00012.

• B.C. Allanach, Tyler Corbett, and Maeve Madigan. Sensitivity of Future Hadron Colliders to Leptoquark Pair Production in the Di-Muon Di-Jets Channel. Eur. Phys. J. C, 80(2):170, 2020. arXiv:1911.04455, doi:10.1140/epjc/s10052-020-7722-3.

• Simone Amoroso, Deepak Kar, and Matthias Schott. How to discover QCD Instantons at the LHC. In Topological Effects in the Standard Model: Instantons, Sphalerons and Beyond at LHC. 12 2020. arXiv:2012.09120.

• Andrei Angelescu, Damir Bečirević, Darius A. Faroughy, Florentin Jaffredo, and Olcyr Sumensari. On the single leptoquark solutions to the $B$-physics anomalies. 3 2021. arXiv:2103.12504.

• Andrei Angelescu, Darius A. Faroughy, and Olcyr Sumensari. Lepton Flavor Violation and Dilepton Tails at the LHC. Eur. Phys. J. C, 80(7):641, 2020. arXiv:2002.05684, doi:10.1140/epjc/s10052-020-8210-5.

• Jack Y. Araz and others. Proceedings of the second MadAnalysis 5 workshop on LHC recasting in Korea. Mod. Phys. Lett. A, 36(01):2102001, 2021. arXiv:2101.02245, doi:10.1142/S0217732321020016.

• Johann Brehmer, Felix Kling, Irina Espejo, and Kyle Cranmer. MadMiner: Machine learning-based inference for particle physics. Comput. Softw. Big Sci., 4(1):3, 2020. arXiv:1907.10621, doi:10.1007/s41781-020-0035-2.

• G. Brooijmans and others. Les Houches 2019 Physics at TeV Colliders: New Physics Working Group Report. In 2020. arXiv:2002.12220.

• Rodolfo Capdevilla, Federico Meloni, Rosa Simoniello, and Jose Zurita. Hunting wino and higgsino dark matter at the muon collider with disappearing tracks. 2 2021. arXiv:2102.11292.

• Vincenzo Cirigliano, Kaori Fuyuto, Christopher Lee, Emanuele Mereghetti, and Bin Yan. Charged Lepton Flavor Violation at the EIC. 2 2021. arXiv:2102.06176.

• Belle II Collaboration. Search for $B^+\to K^+\nu \bar \nu $ decays using an inclusive tagging method at Belle II. 4 2021. arXiv:2104.12624.

• Matthew Feickert, Lukas Heinrich, and Giordon Stark. Likelihood preservation and statistical reproduction of searches for new physics. EPJ Web Conf., 2020. doi:10.1051/epjconf/202024506017.

• Lukas Heinrich, Holger Schulz, Jessica Turner, and Ye-Ling Zhou. Constraining A₄ Leptonic Flavour Model Parameters at Colliders and Beyond. 2018. arXiv:1810.05648.

• Charanjit K. Khosa, Sabine Kraml, Andre Lessa, Philipp Neuhuber, and Wolfgang Waltenberger. SModelS database update v1.2.3. LHEP, 158:2020, 5 2020. arXiv:2005.00555, doi:10.31526/lhep.2020.158.

• Jeffrey Krupa and others. GPU coprocessors as a service for deep learning inference in high energy physics. 7 2020. arXiv:2007.10359.

• Wolfgang Waltenberger, André Lessa, and Sabine Kraml. Artificial Proto-Modelling: Building Precursors of a Next Standard Model from Simplified Model Results. 12 2020. arXiv:2012.12246.

• ATLAS Collaboration. Reproducing searches for new physics with the ATLAS experiment through publication of full statistical likelihoods. Geneva, Aug 2019. URL: https://cds.cern.ch/record/2684863.

• ATLAS Collaboration. Search for new phenomena in events with two opposite-charge leptons, jets and missing transverse momentum in $pp$ collisions at $\sqrt s = 13$ TeV with the ATLAS detector. 2 2021. arXiv:2102.01444.

Published Likelihoods

Updating list of HEPData entries for publications using HistFactory JSON likelihoods:

• Search for bottom-squark pair production with the ATLAS detector in final states containing Higgs bosons, $b$-jets and missing transverse momentum. 2019. URL: https://doi.org/10.17182/hepdata.89408, doi:10.17182/hepdata.89408.

• Search for chargino-neutralino production with mass splittings near the electroweak scale in three-lepton final states in $\sqrt s$ = 13 TeV $pp$ collisions with the ATLAS detector. 2019. URL: https://doi.org/10.17182/hepdata.91127, doi:10.17182/hepdata.91127.

• Search for direct stau production in events with two hadronic τ-leptons in $\sqrt s = 13$ TeV $pp$ collisions with the ATLAS detector. 2019. URL: https://doi.org/10.17182/hepdata.92006, doi:10.17182/hepdata.92006.

• Search for direct production of electroweakinos in final states with one lepton, missing transverse momentum and a Higgs boson decaying into two $b$-jets in (pp) collisions at $\sqrt s=13$ TeV with the ATLAS detector. 2020. URL: https://doi.org/10.17182/hepdata.90607.v2, doi:10.17182/hepdata.90607.v2.

• Search for displaced leptons in $\sqrt s = 13$ TeV $pp$ collisions with the ATLAS detector. 2020. URL: https://doi.org/10.17182/hepdata.98796, doi:10.17182/hepdata.98796.

• Search for squarks and gluinos in final states with same-sign leptons and jets using 139 fb$^-1$ of data collected with the ATLAS detector. 2020. URL: https://doi.org/10.17182/hepdata.91214.v3, doi:10.17182/hepdata.91214.v3.

• Search for trilepton resonances from chargino and neutralino pair production in $\sqrt s$ = 13 TeV $pp$ collisions with the ATLAS detector. 2020. URL: https://doi.org/10.17182/hepdata.99806, doi:10.17182/hepdata.99806.

Roadmap (2019-2020)

This is the pyhf 2019 into 2020 Roadmap (Issue #561).

Overview and Goals

We will follow loosely Seibert’s Heirarchy of Needs

(Stan Seibert, SciPy 2019)

As a general overview that will include:

• Improvements to docs

  • Add lots of examples

  • Add at least 5 well documented case studies

• Issue cleanup

• Adding core feature support

• “pyhf evolution”: integration with columnar data analysis systems

• GPU support and testing

• Publications

  • Submit pyhf to JOSS

  • Submit pyhf to pyOpenSci

  • Start pyhf paper in 2020

• Align with IRIS-HEP Analysis Systems NSF milestones

Time scale

The roadmap will be executed over mostly Quarter 3 of 2019 through Quarter 1 of 2020, with some projects continuing into Quarter 2 of 2020

• 2019-Q3

• 2019-Q4

• 2020-Q1

• (2020-Q2)

Roadmap

1. Documentation and Deployment

   • Add docstrings to all functions and classes (Issues #38, #349) [2019-Q3]

   • Greatly revise and expand examples (Issues #168, #202, #212, #325, #342, #349, #367) [2019-Q3 → 2019-Q4]

     • Add small case studies with published sbottom likelihood from HEPData

   • Move to scikit-hep GitHub organization [2019-Q3]

   • Develop a release schedule/criteria [2019-Q4]

   • Automate deployment with [STRIKEOUT:Azure pipeline (talk with Henry Schreiner) (Issue #517)] GitHub Actions (Issue #508) [2019-Q3]

   • Finalize logo and add it to website (Issue #453) [2019-Q3 → 2019-Q4]

   • Write submission to JOSS (Issue #502) and write submission to pyOpenSci [2019-Q4 → 2020-Q2]

   • Contribute to IRIS-HEP Analysis Systems Milestones “Initial roadmap for ecosystem coherency” and “Initial roadmap for high-level cyberinfrastructure components of analysis system” [2019-Q4 → 2020-Q2]

2. Revision and Maintenance

   • Add tests using HEPData published sbottom likelihoods (Issue #518) [2019-Q3]

   • Add CI with GitHub Actions and Azure Pipelines (PR #527, Issue #517) [2019-Q3]

   • Investigate rewrite of pytest fixtures to use modern pytest (Issue #370) [2019-Q3 → 2019-Q4]

   • Factorize out the statistical fitting portion into pyhf.infer (PR #531) [2019-Q3 → 2019-Q4]

   • Bug squashing at large [2019-Q3 → 2020-Q2]

     • Unexpected use cases (Issues #324, #325, #529)

     • Computational edge cases (Issues #332, #445)

   • Make sure that all backends reproduce sbottom results [2019-Q4 → 2020-Q2]

3. Development

   • Batch support (PR #503) [2019-Q3]

   • Add ParamViewer support (PR #519) [2019-Q3]

   • Add setting of NPs constant/fixed (PR #653) [2019-Q3]

   • Implement pdf as subclass of distributions (PR #551) [2019-Q3]

   • Add sampling with toys (PR #558) [2019-Q3]

   • Make general modeling choices (e.g., Issue #293) [2019-Q4 → 2020-Q1]

   • Add “discovery” test stats (p0) (PR #520) [2019-Q4 → 2020-Q1]

   • Add better Model creation [2019-Q4 → 2020-Q1]

   • Add background model support (Issues #514, #946) [2019-Q4 → 2020-Q1]

   • Develop interface for the optimizers similar to tensor/backend (Issue #754, PR #951) [2019-Q4 → 2020-Q1]

   • Migrate to TensorFlow v2.0 (PR #541) [2019-Q4]

   • Drop Python 2.7 support at end of 2019 (Issue #469) [2019-Q4 (last week of December 2019)]

   • Finalize public API [2020-Q1]

   • Integrate pyfitcore/Statisfactory API [2020-Q1]

4. Research

   • Add pyfitcore/Statisfactory integrations (Issue #344, zfit Issue 120) [2019-Q4]

   • Hardware acceleration scaling studies (Issues #93, #301) [2019-Q4 → 2020-Q1]

   • Speedup through Numba (Issue #364) [2019-Q3 → 2019-Q4]

   • Dask backend (Issue #259) [2019-Q3 → 2020-Q1]

   • Attempt to use pyhf as fitting tool for full Analysis Systems pipeline test in early 2020 [2019-Q4 → 2020-Q1]

   • pyhf should satisfy IRIS-HEP Analysis Systems Milestone “GPU/accelerator-based implementation of statistical and other appropriate components” [2020-Q1 → 2020-Q2] and contributes to “Benchmarking and assessment of prototype analysis system components” [2020-Q3 → 2020-Q4]

Roadmap as Gantt Chart

Presentations During Roadmap Timeline

• Talk at IRIS-HEP Institute Retreat (September 12-13th, 2019)

• Talk at PyHEP 2019 (October 16-18th, 2019)

• Talk at CHEP 2019 (November 4-8th, 2019)

• Poster at CHEP 2019 (November 4-8th, 2019)

Release Notes

v0.6.1

This is a patch release from v0.6.0 → v0.6.1.

Important Notes

• As a result of changes to the default behavior of torch.distributions in PyTorch v1.8.0, accommodating changes have been made in the underlying implementations for pyhf.tensor.pytorch_backend.pytorch_backend(). These changes require a new lower bound of torch v1.8.0 for use of the PyTorch backend.

Fixes

• In the PyTorch backend the validate_args kwarg is used with torch.distributions to ensure a continuous approximation of the Poisson distribution in torch v1.8.0+.

Features

Python API

• The solver_options kwarg can be passed to the pyhf.optimize.opt_scipy.scipy_optimizer() optimizer for additional configuration of the minimization. See scipy.optimize.show_options() for additional options of optimization solvers.

• The torch API is now used to provide the implementations of the ravel, tile, and outer tensorlib methods for the PyTorch backend.

v0.6.0

This is a minor release from v0.5.4 → v0.6.0.

Important Notes

• Please note this release has API breaking changes and carefully read these notes while updating your code to the v0.6.0 API. Perhaps most relevant is the changes to the pyhf.infer.hypotest() API, which now uses a calctype argument to differentiate between using an asymptotic calculator or a toy calculator, and a test_stat kwarg to specify which test statistic the calculator should use, with 'qtilde', corresponding to pyhf.infer.test_statistics.qmu_tilde(), now the default option. It also relies more heavily on using kwargs to pass options through to the optimizer.

• Following the recommendations of NEP 29 — Recommend Python and NumPy version support as a community policy standard pyhf v0.6.0 drops support for Python 3.6. PEP 494 – Python 3.6 Release Schedule also notes that Python 3.6 will be end of life in December 2021, so pyhf is moving forward with a minimum required runtime of Python 3.7.

• Support for the discovery test statistic, \(q_{0}\), has now been added through the pyhf.infer.test_statistics.q0() API.

• Support for pseudoexperiments (toys) has been added through the pyhf.infer.calculators.ToyCalculator() API. Please see the corresponding example notebook for more detailed exploration of the API.

• The minuit extra, python -m pip install pyhf[minuit], now uses and requires the iminuit v2.X release series and API. Note that iminuit v2.X can result in slight differences in minimization results from iminuit v1.X.

• The documentation will now be versioned with releases on ReadTheDocs. Please use pyhf.readthedocs.io to access the documentation for the latest stable release of pyhf.

• pyhf is transtioning away from Stack Overflow to GitHub Discussions for resolving user questions not covered in the documentation. Please check the GitHub Discussions page to search for discussions addressing your questions and to open up a new discussion if your question is not covered.

• pyhf has published a paper in the Journal of Open Source Software. Please make sure to include the paper reference in all citations of pyhf, as documented in the Use and Citations section of the documentation.

Fixes

• Fix bug where all extras triggered warning for installation of the contrib extra.

• float-like values are used in division for pyhf.writexml().

• Model.spec now supports building new models from existing models.

• \(p\)-values are now reported based on their quantiles, instead of interpolating test statistics and converting to \(p\)-values.

• Namespace collisions between uproot3 and uproot/uproot4 have been fixed for the xmlio extra.

• The normsys modifier now uses the pyhf.interpolators.code4 interpolation method by default.

• The histosys modifier now uses the pyhf.interpolators.code4p interpolation method by default.

Features

Python API

• The tensorlib API now supports a tensorlib.to_numpy and tensorlib.ravel API.

• The pyhf.infer.calculators.ToyCalculator() API has been added to support pseudoexperiments (toys).

• The empirical test statistic distribution API has been added to help support the ToyCalculator API.

• Add a tolerance kwarg to the optimizer API to set a float value as a tolerance for termination of the fit.

• The pyhf.optimize.opt_minuit.minuit_optimizer() optimizer now can return correlations of the fitted parameters through use of the return_correlation Boolean kwarg.

• Add the pyhf.utils.citation API to get a str of the preferred BibTeX entry for citation of the version of pyhf installed. See the example for the CLI API for more information.

• The pyhf.infer.hypotest() API now uses a calctype argument to differentiate between using an asymptotic calculator or a toy calculator, and a test_stat kwarg to specify which test statistic to use. It also relies more heavily on using kwargs to pass options through to the optimizer.

• The default test_stat kwarg for pyhf.infer.hypotest() and the calculator APIs is 'qtilde', which corresponds to the alternative test statistic pyhf.infer.test_statistics.qmu_tilde().

• The return type of \(p\)-value like functions is now a 0-dimensional tensor (with shape ()) instead of a float. This is required to support end-to-end automatic differentiation in future releases.

CLI API

• The CLI API now supports a --citation or --cite option to print the preferred BibTeX entry for citation of the version of pyhf installed.

$ pyhf --citation
@software{pyhf,
  author = {Lukas Heinrich and Matthew Feickert and Giordon Stark},
  title = "{pyhf: v0.6.0}",
  version = {0.6.0},
  doi = {10.5281/zenodo.1169739},
  url = {https://github.com/scikit-hep/pyhf},
}

@article{pyhf_joss,
  doi = {10.21105/joss.02823},
  url = {https://doi.org/10.21105/joss.02823},
  year = {2021},
  publisher = {The Open Journal},
  volume = {6},
  number = {58},
  pages = {2823},
  author = {Lukas Heinrich and Matthew Feickert and Giordon Stark and Kyle Cranmer},
  title = {pyhf: pure-Python implementation of HistFactory statistical models},
  journal = {Journal of Open Source Software}
}

Contributors

v0.6.0 benefited from contributions from:

• Alexander Held

• Marco Gorelli

• Pradyumna Rahul K

• Eric Schanet

• Henry Schreiner

v0.5.4

This is a patch release from v0.5.3 → v0.5.4.

Fixes

• Require uproot3 instead of uproot v3.X releases to avoid conflicts when uproot4 is installed in an environment with uproot v3.X installed and namespace conflicts with uproot-methods. Adoption of uproot3 in v0.5.4 will ensure v0.5.4 works far into the future if XML and ROOT I/O through uproot is required.

Example:

Without the v0.5.4 patch release there is a regression in using uproot v3.X and uproot4 in the same environment (which was swiftly identified and patched by the fantastic uproot team)

$ python -m pip install "pyhf[xmlio]<0.5.4"
$ python -m pip list | grep "pyhf\|uproot"
pyhf           0.5.3
uproot         3.13.1
uproot-methods 0.8.0
$ python -m pip install uproot4
$ python -m pip list | grep "pyhf\|uproot"
pyhf           0.5.3
uproot         4.0.0
uproot-methods 0.8.0
uproot4        4.0.0

this is resolved in v0.5.4 with the requirement of uproot3

$ python -m pip install "pyhf[xmlio]>=0.5.4"
$ python -m pip list | grep "pyhf\|uproot"
pyhf            0.5.4
uproot3         3.14.1
uproot3-methods 0.10.0
$ python -m pip install uproot4 # or uproot
$ python -m pip list | grep "pyhf\|uproot"
pyhf            0.5.4
uproot          4.0.0
uproot3         3.14.1
uproot3-methods 0.10.0
uproot4         4.0.0

v0.5.3

This is a patch release from v0.5.2 → v0.5.3.

Fixes

• Workspaces are now immutable

• ShapeFactor support added to XML reading and writing

• An error is raised if a fit initialization parameter is outside of its bounds (preventing hypotest with POI outside of bounds)

Features

Python API

• Inverting hypothesis tests to get upper limits now has an API with pyhf.infer.intervals.upperlimit

• Building workspaces from a model and data added with pyhf.workspace.build

CLI API

• Added CLI API for pyhf.infer.fit: pyhf fit

• pyhf combine now allows for merging channels: pyhf combine --merge-channels --join <join option>

• Added utility to download archived pyhf pallets (workspaces + patchsets) to contrib module: pyhf contrib download

Contributors

v0.5.3 benefited from contributions from:

• Karthikeyan Singaravelan

Contributors

pyhf is openly developed and benefits from the contributions and feedback from its users. The pyhf dev team would like to thank all contributors to the project for their support and help. Thank you!

Contributors include:

• Jessica Forde

• Ruggero Turra

• Tadej Novak

• Frank Sauerburger

• Lars Nielsen

• Kanishk Kalra

• Nikolai Hartmann

• Alexander Held

• Karthikeyan Singaravelan

• Marco Gorelli

• Pradyumna Rahul K

• Eric Schanet

• Henry Schreiner

View me on GitHub

Warning: This is a development version. The latest stable version is at ReadTheDocs.

pure-python fitting/limit-setting/interval estimation HistFactory-style

The HistFactory p.d.f. template [CERN-OPEN-2012-016] is per-se independent of its implementation in ROOT and sometimes, it’s useful to be able to run statistical analysis outside of ROOT, RooFit, RooStats framework.

This repo is a pure-python implementation of that statistical model for multi-bin histogram-based analysis and its interval estimation is based on the asymptotic formulas of “Asymptotic formulae for likelihood-based tests of new physics” [arXiv:1007.1727]. The aim is also to support modern computational graph libraries such as PyTorch and TensorFlow in order to make use of features such as autodifferentiation and GPU acceleration.

Hello World

This is how you use the pyhf Python API to build a statistical model and run basic inference:

>>> import pyhf
>>> model = pyhf.simplemodels.hepdata_like(
...     signal_data=[12.0, 11.0], bkg_data=[50.0, 52.0], bkg_uncerts=[3.0, 7.0]
... )
>>> data = [51, 48] + model.config.auxdata
>>> test_mu = 1.0
>>> CLs_obs, CLs_exp = pyhf.infer.hypotest(
...     test_mu, data, model, test_stat="qtilde", return_expected=True
... )
>>> print(f"Observed: {CLs_obs}, Expected: {CLs_exp}")
Observed: 0.05251497423736956, Expected: 0.06445320535890459

Alternatively the statistical model and observational data can be read from its serialized JSON representation (see next section).

>>> import pyhf
>>> import requests
>>> wspace = pyhf.Workspace(requests.get("https://git.io/JJYDE").json())
>>> model = wspace.model()
>>> data = wspace.data(model)
>>> test_mu = 1.0
>>> CLs_obs, CLs_exp = pyhf.infer.hypotest(
...     test_mu, data, model, test_stat="qtilde", return_expected=True
... )
>>> print(f"Observed: {CLs_obs}, Expected: {CLs_exp}")
Observed: 0.3599840922126626, Expected: 0.3599840922126626

Finally, you can also use the command line interface that pyhf provides

$ cat << EOF  | tee likelihood.json | pyhf cls
{
    "channels": [
        { "name": "singlechannel",
          "samples": [
            { "name": "signal",
              "data": [12.0, 11.0],
              "modifiers": [ { "name": "mu", "type": "normfactor", "data": null} ]
            },
            { "name": "background",
              "data": [50.0, 52.0],
              "modifiers": [ {"name": "uncorr_bkguncrt", "type": "shapesys", "data": [3.0, 7.0]} ]
            }
          ]
        }
    ],
    "observations": [
        { "name": "singlechannel", "data": [51.0, 48.0] }
    ],
    "measurements": [
        { "name": "Measurement", "config": {"poi": "mu", "parameters": []} }
    ],
    "version": "1.0.0"
}
EOF

which should produce the following JSON output:

{
   "CLs_exp": [
      0.0026062609501074576,
      0.01382005356161206,
      0.06445320535890459,
      0.23525643861460702,
      0.573036205919389
   ],
   "CLs_obs": 0.05251497423736956
}

What does it support

Implemented variations:

• ☑ HistoSys

• ☑ OverallSys

• ☑ ShapeSys

• ☑ NormFactor

• ☑ Multiple Channels

• ☑ Import from XML + ROOT via uproot

• ☑ ShapeFactor

• ☑ StatError

• ☑ Lumi Uncertainty

• ☑ Non-asymptotic calculators

Computational Backends:

• ☑ NumPy

• ☑ PyTorch

• ☑ TensorFlow

• ☑ JAX

Optimizers:

• ☑ SciPy (scipy.optimize)

• ☑ MINUIT (iminuit)

All backends can be used in combination with all optimizers. Custom user backends and optimizers can be used as well.

Todo

• ☐ StatConfig

results obtained from this package are validated against output computed from HistFactory workspaces

A one bin example

import pyhf
import numpy as np
import matplotlib.pyplot as plt
from pyhf.contrib.viz import brazil

pyhf.set_backend("numpy")
model = pyhf.simplemodels.hepdata_like(
    signal_data=[10.0], bkg_data=[50.0], bkg_uncerts=[7.0]
)
data = [55.0] + model.config.auxdata

poi_vals = np.linspace(0, 5, 41)
results = [
    pyhf.infer.hypotest(
        test_poi, data, model, test_stat="qtilde", return_expected_set=True
    )
    for test_poi in poi_vals
]

fig, ax = plt.subplots()
fig.set_size_inches(7, 5)
brazil.plot_results(poi_vals, results, ax=ax)
fig.show()

pyhf

ROOT

A two bin example

import pyhf
import numpy as np
import matplotlib.pyplot as plt
from pyhf.contrib.viz import brazil

pyhf.set_backend("numpy")
model = pyhf.simplemodels.hepdata_like(
    signal_data=[30.0, 45.0], bkg_data=[100.0, 150.0], bkg_uncerts=[15.0, 20.0]
)
data = [100.0, 145.0] + model.config.auxdata

poi_vals = np.linspace(0, 5, 41)
results = [
    pyhf.infer.hypotest(
        test_poi, data, model, test_stat="qtilde", return_expected_set=True
    )
    for test_poi in poi_vals
]

fig, ax = plt.subplots()
fig.set_size_inches(7, 5)
brazil.plot_results(poi_vals, results, ax=ax)
fig.show()

pyhf

ROOT

Installation

To install pyhf from PyPI with the NumPy backend run

python -m pip install pyhf

and to install pyhf with all additional backends run

python -m pip install pyhf[backends]

or a subset of the options.

To uninstall run

python -m pip uninstall pyhf

Questions

If you have a question about the use of pyhf not covered in the documentation, please ask a question on the GitHub Discussions.

If you believe you have found a bug in pyhf, please report it in the GitHub Issues. If you’re interested in getting updates from the pyhf dev team and release announcements you can join the pyhf-announcements mailing list.

Citation

As noted in Use and Citations, the preferred BibTeX entry for citation of pyhf includes both the Zenodo archive and the JOSS paper:

@software{pyhf,
  author = {Lukas Heinrich and Matthew Feickert and Giordon Stark},
  title = "{pyhf: v0.6.1}",
  version = {0.6.1},
  doi = {10.5281/zenodo.1169739},
  url = {https://github.com/scikit-hep/pyhf},
}

@article{pyhf_joss,
  doi = {10.21105/joss.02823},
  url = {https://doi.org/10.21105/joss.02823},
  year = {2021},
  publisher = {The Open Journal},
  volume = {6},
  number = {58},
  pages = {2823},
  author = {Lukas Heinrich and Matthew Feickert and Giordon Stark and Kyle Cranmer},
  title = {pyhf: pure-Python implementation of HistFactory statistical models},
  journal = {Journal of Open Source Software}
}

Authors

pyhf is openly developed by Lukas Heinrich, Matthew Feickert, and Giordon Stark.

Please check the contribution statistics for a list of contributors.

Milestones

• 2020-07-28: 1000 GitHub issues and pull requests. (See PR #1000)

Acknowledgements

Matthew Feickert has received support to work on pyhf provided by NSF cooperative agreement OAC-1836650 (IRIS-HEP) and grant OAC-1450377 (DIANA/HEP).

Indices and tables

• Index

• Module Index

• Search Page