statistics package

Submodules

statistics.analysesCombinations module

class statistics.analysesCombinations.AnaCombLikelihoodComputer(theoryPredictions: list, deltas_rel=None)

Bases: object

constructor.

Parameters

• theoryPredictions – the List of theory predictions

• deltas_rel – relative uncertainty in signal (float). Default value is 20%.

CLs(mu: float = 1.0, evaluationType: NllEvalType = NllEvalType.observed, return_type: str = 'CLs', nSigma: int = 0) → float

Compute the exclusion confidence level of the model

Parameters

• mu – compute for the parameter of interest mu

• evaluationType – one of: observed, apriori, aposteriori

• return_type – (Text) can be “CLs-alpha”, “1-CLs”, “CLs” CLs-alpha: returns CLs - 0.05 alpha-CLs: returns 0.05 - CLs 1-CLs: returns 1-CLs value CLs: returns CLs value

getCLsRootFunc(evaluationType: NllEvalType = NllEvalType.observed, allowNegativeSignals: bool = False, nSigma: int = 0) → Tuple[float, float, Callable]

Obtain the function “CLs-alpha[0.05]” whose root defines the upper limit, plus mu_hat and sigma_mu :param nSigma: the upper limit for central value (0), + 1 sigma, - 1 sigma, etc. For error bands. :param evaluationType: one of: observed, apriori, aposteriori

getLlhds(muvals, evaluationType: NllEvalType = NllEvalType.observed, normalize: bool = True)

Compute the likelihoods for the individual analyses and the combined likelihood. Returns a dictionary with the analysis IDs as keys and the likelihood values as values.

Parameters

muvals – List with values for the signal strenth for which the

likelihoods must be evaluated. :param evaluationType: one of: observed, apriori, aposteriori :param normalize: If True normalizes the likelihood by its integral over muvals.

getUpperLimitOnMu(evaluationType: NllEvalType = NllEvalType.observed, allowNegativeSignals=False, nSigma: int = 0) → float

get upper limit on signal strength multiplier, i.e. value for mu for which CLs = 0.95

Parameters

• evaluationType – one of: observed, apriori, aposteriori

• nSigma – the upper limit for central value (0),

• 1 sigma, - 1 sigma, etc. For error bands.

Returns

upper limit on signal strength multiplier mu

getUpperLimitOnSigmaTimesEff(evaluationType: NllEvalType = NllEvalType.observed, allowNegativeSignals: bool = False, nSigma: int = 0) → Unum

upper limit on the fiducial cross section sigma times efficiency,

summed over all signal regions, i.e. sum_i xsec^prod_i eff_i obtained from the defined Data (using the signal prediction for each signal region/dataset), by using the q_mu test statistic from the CCGV paper (arXiv:1007.1727).

Parameters

• evaluationType – one of: observed, apriori, aposteriori

• nSigma – the upper limit for central value (0),

• 1 sigma, - 1 sigma, etc. For error bands.

Returns

upper limit on fiducial cross section

likelihood(mu: float = 1.0, evaluationType: NllEvalType = NllEvalType.observed, return_nll: bool = False, asimov: Union[None, float] = None) → float

Compute the likelihood at a given mu

Parameters

• mu – signal strength

• evaluationType – one of: observed, apriori, aposteriori

• return_nll – if True, return negative log likelihood, else likelihood

• asimov – if not None, compute llhd for asimov data with mu=asimov

lmax(allowNegativeSignals: bool = False, evaluationType: NllEvalType = NllEvalType.observed, return_nll: bool = False, asimov: Union[None, float] = None) → Optional[Dict]

find muhat and lmax.

Parameters

• allowNegativeSignals – if true, then also allow for negative values

• evaluationType – one of: observed, apriori, aposteriori

• return_nll – if true, return negative log max likelihood instead of lmax

Returns

mu_hat, i.e. the maximum likelihood estimate of mu, if extended output is requested, it returns a dictionary with mu_hat, sigma_mu – the standard deviation around mu_hat, and lmax, i.e. the likelihood at mu_hat

nll(mu: float = 1.0, evaluationType: NllEvalType = NllEvalType.observed, asimov: Union[None, float] = None) → float

Compute the negative log likelihood at a given mu. For now this is just a proxy method to ease transition of codes that depend on SModelS.

Parameters

• mu – signal strength

• evaluationType – one of: observed, apriori, aposteriori

• asimov – if not None, compute llhd for asimov data with mu=asimov

statistics.basicStats module

statistics.basicStats.CLsfromNLL(nllA: float, nll0A: float, nll: float, nll0: float, big_muhat: bool, return_type: str = 'CLs-alpha', nSigma: int = 0) → float

compute CLs (or similar) from the NLLs :param nllA: negative log likelihood for Asimov data :param nll0A: negative log likelihood at muhat for Asimov data :param nll: negative log likelihood :param nll0: negative log likelihood at muhat :param big_muhat: true if muhat>mu :param return_type: (Text) one of “CLs-alpha”, “alpha-CLs”, “1-CLs”, or “CLs”: ========== ======================= CLs-alpha: returns CLs - 0.05 alpha-CLs: returns 0.05 - CLs 1-CLs: returns 1-CLs value CLs: returns “raw” CLs value ========== ======================= :param nSigma: compute CLs not for central value but for this number of sigmas (positive or negative), see Equations 86 - 89 in the CCGV paper.

Returns

Cls-type value, see above

class statistics.basicStats.NllEvalType(value, names=None, *, module=None, qualname=None, type=None, start=1, boundary=None)

Bases: Enum

an enum to account for the different types of likelihood values: observed, a priori evaluationType, a posteriori evaluationType

aposteriori = 1

apriori = 2

classmethod init(evaluationType: Union[str, bool])

get evaluationtype either from a string (e.g. ‘posteriori’) or a bool (true is priori, false is observed)

observed = 0

statistics.basicStats.chi2FromLmax(llhd, lmax)

compute the chi2 from likelihood and lmax

statistics.basicStats.determineBrentBracket(mu_hat, sigma_mu, rootfinder, allowNegative=True)

find a, b for brent bracketing

Parameters

• mu_hat – mu that maximizes likelihood

• sigm_mu – error on mu_hat (not too reliable)

• rootfinder – function that finds the root (usually root_func)

• allowNegative – if False, then do not allow a or b to become negative

Returns

the interval a,b

statistics.exceptions module

exception statistics.exceptions.SModelSStatisticsError(value=None)

Bases: Exception

Class to define SModelS specific errors

statistics.pyhfInterface module

class statistics.pyhfInterface.PyhfData(nsignals: Dict[str, Dict], inputJsons, jsonFiles=None, includeCRs=False, signalUncertainty=None, globalInfo=None)

Bases: object

Holds data for use in pyhf :ivar nsignals: signal predictions dictionary of dictionaries, one for each json file, one entry per signal region :ivar inputJsons: list of json instances :ivar jsonFiles: optional list of json files :ivar nWS: number of workspaces = number of json files :ivar globalInfo: None, or link to globalInfo (for debugging)

checkConsistency()

Check various inconsistencies of the PyhfData attributes

Parameters

zeroSignalsFlag – boolean identifying if all SRs of a single json are empty

createPatchForRegion(region, i_ch, ch, jsName)

getTotalYield()

the total yield in all signal regions

getWSInfo()

Getting informations from the json files

Variables

channelsInfo – list of dictionaries (one dictionary for each json file) containing useful information about the json files - :key signalRegions: list of dictonaries with ‘json path’ and ‘size’ (number of bins) of the ‘signal regions’ channels in the json files - :key otherRegions: list of dictionnaries indicating the path and the name of the control and/or validation region channels

updatePyhfNames(jsName: str, observations: List)

if no pyhf names are given, get them from the ws, in the order of the ws

class statistics.pyhfInterface.PyhfUpperLimitComputer(data, cl=0.95, lumi=None)

Bases: object

Class that computes the upper limit using the jsons files and signal informations in the ‘data’ instance of ‘PyhfData’

Parameters

• data – instance of ‘PyhfData’ holding the signals information

• cl – confdence level at which the upper limit is desired to be computed

Variables

• data – created from data

• nsignals – signal predictions list divided into sublists, one for each json file

• inputJsons – list of input json files as python json instances

• channelsInfo – list of channels information for the json files

• zeroSignalsFlag – list boolean flags in case all signals are zero for a specific json

• nWS – number of workspaces = number of json files

• patches – list of patches to be applied to the inputJsons as python dictionary instances

• workspaces – list of workspaces resulting from the patched inputJsons

• workspaces_expected – list of patched workspaces with observation yields replaced by the expected eones

• cl – created from cl

• scale – scale that is applied to the signal predictions, dynamically changes throughout the upper limit calculation

• alreadyBeenThere – boolean flag that identifies when nsignals accidentally passes twice at two identical values

CLs(mu: float, evaluationType: NllEvalType, return_type: str = 'CLs', workspace_index: Optional[int] = None, nSigma: int = 0) → float

This is the CLs method that heeds self.scale, i.e. you give it the _absolute_ mu

Parameters

• mu – compute for the parameter of interest mu

• evaluationType – if observed, compute observed, apriori: compute a priori expected

• return_type – (Text) can be one of:

“CLs-alpha”, “1-CLs”, “CLs” “alpha-CLs” CLs-alpha: returns CLs - 0.05 alpha-CLs: return 0.05 - CLs 1-CLs: returns 1-CLs value CLs: returns CLs value

backup()

changeChannelName(srInfo)

changes the channel names in the json to match the SModelS name. FIXME this is only a hack for now, should be replaced by a self-respecting dataIdMap in the database itself, that maps the smodels names to the pyhf names explicitly. This method will then go away!

checkPyhfVersion()

check the pyhf version, currently we need 0.6.1+

compute_invhess(x, data, model, index, epsilon=1e-05)

if inv_hess is not given by the optimiser, calculate numerically by evaluating second order partial derivatives using 2 point central finite differences method :param x: parameter values given to pyhf.infer.mle.twice_nll taken from pyhf.infer.mle.fit - optimizer.x (best_fit parameter values) :param data: workspace.data(model) passed to pyhf.infer.mle.fit :param model: model passed to pyhf.infer.mle.fit :param index: index of the POI Note : If len(x) <=5, compute the entire hessian matrix and ind its inverse. Else, compute the hessian at the index of the POI and return its inverse (diagonal approximation) returns the inverse hessian at the index of the poi

exponentiateNLL(twice_nll, doIt)

if doIt, then compute likelihood from nll, else return nll

generateAsimovData(mu: Union[None, float] = 0.0, workspace_index: Optional[int] = None, evaluationType: NllEvalType = NllEvalType.observed) → list

generate asimov data for the model for signal strength mu :param mu: if None, then actually return the original data :returns: tuple with workspace and list with asimov data

getBestCombinationIndex()

find the index of the best evaluationType combination

getUpperLimitOnMu(evaluationType: NllEvalType = NllEvalType.observed, workspace_index=None, nSigma: int = 0)

Compute the upper limit on the signal strength modifier with:

• by default, the combination of the workspaces contained into self.workspaces

• if workspace_index is specified, self.workspace[workspace_index] (useful for computation of the best upper limit)

Parameters

• evaluationType –

  • if set to apriori: uses evaluationType SM backgrounds as signals

  • else: uses ‘self.nsignals’

• workspace_index –

  • if different from ‘None’: index of the workspace to use for upper limit

  • else: choose best combo

Returns

the upper limit at ‘self.cl’ level (0.95 by default)

getUpperLimitOnSigmaTimesEff(evaluationType: NllEvalType = NllEvalType.observed, workspace_index=None, nSigma: int = 0)

Compute the upper limit on the fiducial cross section sigma times efficiency:

• by default, the combination of the workspaces contained into self.workspaces

• if workspace_index is specified, self.workspace[workspace_index] (useful for computation of the best upper limit)

Parameters

• evaluationType –

  • if set to apriori: uses evaluationType SM backgrounds as signals

  • else: uses ‘self.nsignals’

• workspace_index –

  • if different from ‘None’: index of the workspace to use for upper limit

  • else: choose best combo

Returns

the upper limit on sigma times eff at ‘self.cl’ level (0.95 by default)

get_position(name, model)

Parameters

• name – name of the parameter one wants to increase

• model – the pyhf model

Returns

the position of the parameter that has to be modified in order to turn positive the negative total yield

likelihood(mu: float = 1.0, workspace_index: Union[None, int] = None, return_nll: bool = False, evaluationType: NllEvalType = NllEvalType.observed, asimov: Union[None, float] = None)

Returns the value of the likelihood. Inspired by the ‘pyhf.infer.mle’ module but for non-log likelihood

Parameters

• workspace_index – supply index of workspace to use. If None, choose index of best combo

• return_nll – if true, return nll, not llhd

• evaluationType – one of: observed, apriori, aposteriori

lmax(workspace_index=None, return_nll=False, evaluationType: NllEvalType = NllEvalType.observed, allowNegativeSignals=False)

Returns the negative log max likelihood

Parameters

• return_nll – if true, return nll, not llhd

• workspace_index – supply index of workspace to use. If None, choose index of best combo

• evaluationType – one of: observed, apriori, aposteriori

• allowNegativeSignals – if False, then negative nsigs are replaced with 0.

patchMaker()

Method that creates the list of patches to be applied to the self.inputJsons workspaces, one for each region given the self.nsignals and the informations available in self.channelsInfo and the content of the self.inputJsons NB: It seems we need to include the change of the “modifiers” in the patches as well

Returns

the list of patches, one for each workspace

rescale(factor)

Rescales the signal predictions (self.nsignals) and processes again the patches and workspaces

Returns

updated list of patches and workspaces (self.patches, self.workspaces and self.workspaces_expected)

rescaleBgYields(init_pars, workspace, model)

Increase initial value of nuisance parameters until the starting value of the total yield (mu*signal + background) is positive

Parameters

• init_pars – list of initial parameters values one wants to increase in order to turn positive the negative total yields

• workspace – the pyhf workspace

• model – the pyhf model

Returns

the list of initial parameters values that gives positive total yields

restore()

updateWorkspace(workspace_index=None, evaluationType: NllEvalType = NllEvalType.observed)

Small method used to return the appropriate workspace

Parameters

• workspace_index – the index of the workspace to retrieve from the corresponding list

• evaluationType – if False, retuns the unmodified (but patched) workspace. Used for computing observed or aposteriori evaluationType limits. if True, retuns the modified (and patched) workspace, where obs = sum(bkg). Used for computing apriori evaluationType limit.

welcome()

greet the world

wsMaker(apriori=False)

Apply each region patch (self.patches) to his associated json (self.inputJsons) to obtain the complete workspaces

Parameters

apriori –

• If set to ‘True’: Replace the observation data entries of each workspace by the corresponding sum of the evaluationType yields - Else: The observed yields put in the workspace are the ones written in the corresponfing json dictionary

Returns

the list of patched workspaces

statistics.pyhfInterface.guessPyhfType(name: str) → str

given the pyhf analysis region name, guess the type. awkward, phase this out!

statistics.pyhfInterface.setBackend(backend: str) → bool

try to setup backend to <backend>

Parameters

backend – one of: numpy (default), pytorch, jax, tensorflow

Returns

True, if worked, False if failed

statistics.simplifiedLikelihoods module

class statistics.simplifiedLikelihoods.Data(observed, backgrounds, covariance, third_moment=None, nsignal=None, name: str = 'model', deltas_rel: float = 0.2, lumi=None, asimov: Union[None, float] = None)

Bases: object

A very simple observed container to collect all the data needed to fully define a specific statistical model

Parameters

• observed – number of observed events per dataset

• backgrounds – evaluationType bg per dataset

• covariance – uncertainty in background, as a covariance matrix

• nsignal – number of signal events in each dataset

• name – give the model a name, just for convenience

• deltas_rel – the assumed relative error on the signal hypotheses. The default is 20%.

• lumi – luminosity of dataset in 1/fb, or None

• asimov – None if not asimov data, else signal strength of asimov data

convert(obj)

Convert object to numpy arrays. If object is a float or int, it is converted to a one element array.

correlations()

Correlation matrix, computed from covariance matrix. Convenience function.

diagCov()

Diagonal elements of covariance matrix. Convenience function.

generateAsimovData(theta_hat: list, mu: float = 0.0)

generate a model with Asimov data out of the current model :param theta_hat: the profiled nuisance parameters :returns: Data object

isLinear()

Statistical model is linear, i.e. no quadratic term in poissonians

isScalar(obj)

Determine if obj is a scalar (float or int)

nsignals(mu)

Returns the number of evaluationType signal events, for all datasets, given total signal strength mu.

Parameters

mu – Total number of signal events summed over all datasets.

rel_signals(mu)

Returns the number of evaluationType relative signal events, for all datasets, given total signal strength mu. For mu=1, the sum of the numbers = 1.

Parameters

mu – Total number of signal events summed over all datasets.

sandwich()

Sandwich product

totalCovariance(nsig: list)

get the total covariance matrix, taking into account also signal uncertainty for the signal hypothesis <nsig>. If nsig is None, the predefined signal hypothesis is taken.

var_s(nsig: Union[None, list] = None)

The signal variances. Convenience function.

Parameters

nsig – If None, it will use the model evaluationType number of signal events, otherwise will return the variances for the input value using the relative signal uncertainty defined for the model.

zeroSignal()

Is the total number of signal events zero?

class statistics.simplifiedLikelihoods.LikelihoodComputer(data)

Bases: object

Parameters

data – a Data object.

chi2()

Computes the chi2 for a given number of observed events nobs given the predicted background nb, error on this background deltab, evaluationType number of signal events nsig and the relative error on signal (deltas_rel).

Parameters

nsig – number of signal events

Returns

chi2 (float)

d2NLLdMu2(mu, theta_hat, allowZeroHessian=True)

the hessian of the likelihood of mu, at mu, which is the Fisher information which is approximately the inverse of the covariance

Parameters

allowZeroHessian – if false and sum(observed)==0, then replace observed with expected

d2NLLdTheta2(theta)

the Hessian of nll as a function of the thetas. Makes it easier to find the maximum likelihood.

dNLLdMu(mu, theta_hat=None)

d (- ln L)/d mu, if L is the likelihood. The function whose root gives us muhat, i.e. the mu that maximizes the likelihood.

Parameters

• mu – total number of signal events

• theta_hat – array with nuisance parameters, if None then compute them

dNLLdTheta(theta)

the derivative of nll as a function of the thetas. Makes it easier to find the maximum likelihood.

debug_mode = False

extendedOutput(extended_output: bool, default: Union[None, float] = None, return_nll: bool = False)

findAvgr(theta_hat)

from the difference observed - background, find got inital values for lower and upper

findMuHat(allowNegativeSignals=False, extended_output=False, return_nll=False)

Find the most likely signal strength mu via gradient descent given the relative signal strengths in each dataset (signal region).

Parameters

• allowNegativeSignals – if true, then also allow for negative values

• extended_output – if true, return also sigma_mu, the estimate of the error of mu_hat, and lmax, the likelihood at mu_hat

• return_nll – if true, return nll instead of lmax in the extended output

Returns

mu_hat, i.e. the maximum likelihood estimate of mu, if extended output is requested, it returns mu_hat, sigma_mu – the standard deviation around mu_hat, and llhd, the likelihood at mu_hat

findMuHatViaBracketing(allowNegativeSignals=False, extended_output=False, return_nll=False)

Find the most likely signal strength mu via a brent bracketing technique given the relative signal strengths in each dataset (signal region).

Parameters

• allowNegativeSignals – if true, then also allow for negative values

• extended_output – if true, return also sigma_mu, the estimate of the error of mu_hat, and lmax, the likelihood at mu_hat

• return_nll – if true, return nll instead of lmax in the

extended output

Returns

mu_hat, i.e. the maximum likelihood estimate of mu, if extended output is requested, it returns a dictionary with mu_hat, sigma_mu – the standard deviation around mu_hat, and lmax, i.e. the likelihood at mu_hat

findThetaHat(mu: float)

Compute nuisance parameters theta that maximize our likelihood (poisson*gauss).

generateAsimovComputer(mu: float = 0.0)

given the current computer with its model, generate a model with asimov data for a given mu. :param mu: signal strength to create Asimov data for.

Returns

LikelihoodComputer, with theta_hat as data member

getSigmaMu(mu, theta_hat)

Get an estimate for the standard deviation of mu at <mu>, from the inverse hessian

getThetaHat(nobs, nb, mu, covb, max_iterations)

Compute nuisance parameter theta that maximizes our likelihood (poisson*gauss) – by setting dNLL/dTheta to zero

Parameters

mu – signal strength

Returns

theta_hat

likelihood(mu: float, return_nll: bool = False, evaluationType: NllEvalType = NllEvalType.observed, asimov: Union[None, float] = None)

compute the profiled likelihood for mu.

Parameters

• mu – float Parameter of interest, signal strength

• return_nll – if true, return nll instead of likelihood

Returns

profile likelihood and error code (0=no error)

llhdOfTheta(theta, nll=True)

likelihood for nuicance parameters theta, given signal strength self.mu. notice, by default it returns nll

Parameters

• theta – nuisance parameters

• nll – if True, compute negative log likelihood

lmax(return_nll=False, allowNegativeSignals=False)

convenience function, computes likelihood for nsig = nobs-nbg,

Parameters

• return_nll – return nll instead of likelihood

• allowNegativeSignals – if False, then negative nsigs are replaced with 0.

transform(evaluationType: NllEvalType)

replace the actual observations with backgrounds, if evaluationType is True or “posteriori”

class statistics.simplifiedLikelihoods.UpperLimitComputer(likelihoodComputer, cl: float = 0.95)

Bases: object

Parameters

cl – desired quantile for limits

CLs(mu: float = 1.0, evaluationType: NllEvalType = NllEvalType.observed, trylasttime: bool = False, return_type: str = 'CLs') → float

Compute the exclusion confidence level of the model (1-CLs).

Parameters

• model – statistical model

• evaluationType – one of: observed, apriori, aposteriori

• trylasttime – if True, then dont try extra

• return_type – (Text) can be “CLs-alpha”, “1-CLs”, “CLs” CLs-alpha: returns CLs - 0.05 (alpha) 1-CLs: returns 1-CLs value CLs: returns CLs value

debug_mode = False

getCLsRootFunc(evaluationType: NllEvalType = NllEvalType.observed, trylasttime: Optional[bool] = False, nSigma: int = 0) → Tuple

Obtain the function “CLs-alpha[0.05]” whose root defines the upper limit, plus mu_hat and sigma_mu

Parameters

• model – statistical model

• evaluationType – one of: observed, apriori, aposteriori

• trylasttime – if True, then dont try extra

• nSigma – the upper limit for central value (0),

• 1 sigma, - 1 sigma, etc.

For error bands. :return: mu_hat, sigma_mu, CLs-alpha

getUpperLimitOnMu(evaluationType: NllEvalType = NllEvalType.observed, trylasttime: bool = False, nSigma: int = 0) → Union[None, float]

upper limit on the signal strength multiplier mu obtained from the defined Data (using the signal prediction for each signal regio/dataset), by using the q_mu test statistic from the CCGV paper (arXiv:1007.1727).

Parameters

• evaluationType – one of: observed, apriori, aposteriori.

• trylasttime – if True, then dont try extra.

• nSigma – the upper limit for central value (0),

• 1 sigma, - 1 sigma, etc.

For error bands. :returns: upper limit on the signal strength multiplier mu

getUpperLimitOnSigmaTimesEff(evaluationType: NllEvalType = NllEvalType.observed, trylasttime: bool = False, nSigma: int = 0) → Unum

upper limit on the fiducial cross section sigma times efficiency,

summed over all signal regions, i.e. sum_i xsec^prod_i eff_i obtained from the defined Data (using the signal prediction for each signal region/dataset), by using the q_mu test statistic from the CCGV paper (arXiv:1007.1727).

Parameters

• evaluationType – one of: observed, apriori, aposteriori

• trylasttime – if True, then dont try extra

• nSigma – the upper limit for central value (0),

• 1 sigma, - 1 sigma, etc. For error bands.

Returns

upper limit on fiducial cross section

statistics.speyPyhf module

class statistics.speyPyhf.SpeyPyhfData(nsignals: list, inputJsons, jsonFiles=None, includeCRs: bool = False)

Bases: object

Holds data for use in pyhf :ivar nsignals: signal predictions list divided into sublists, one for each

json file

Variables

• inputJsons – list of json instances

• jsonFiles – optional list of json files

• nWS – number of workspaces = number of json files

channelsInfo

checkConsistency()

Check various inconsistencies of the PyhfData attributes

Parameters

zeroSignalsFlag – boolean identifying if all SRs of a single json are empty

classmethod createDataObject(dataset, nsig: list)

an object creator method

errorFlag

getTotalYield()

the total yield in all signal regions

getWSInfo()

Getting information from the json files

Variables

channelsInfo – list of dictionaries (one dictionary for each json

file) containing useful information about the json files

• key signalRegions

  list of dictonaries with ‘json path’ and

  ‘size’ (number of bins) of the ‘signal regions’ channels in the json files

• key otherRegions

  list of strings indicating the path to the

  control and validation region channels

includeCRs

inputJsons

jsonFiles

nWS

nsignals

patchMaker()

Method that creates the list of patches to be applied to the self.inputJsons workspaces, one for each region given the self.nsignals and the information available in self.channelsInfo and the content of the self.inputJsons NB: It seems we need to include the change of the “modifiers” in the patches as well

Returns

the list of patches, one for each workspace

totalYield

wsMaker(apriori=False)

Apply each region patch (self.patches) to his associated json

(self.inputJsons) to obtain the complete workspaces

Parameters

apriori –

• If set to True: Replace the observation data

entries of each workspace by the corresponding sum of the expected yields - Else: The observed yields put in the workspace are the ones written in the corresponfing json dictionary

Returns

the list of patched workspaces

zeroSignalsFlag

statistics.speyTools module

statistics.statsTools module

class statistics.statsTools.StatsComputer(dataObject: Union[DataSet, CombinedDataSet, list], dataType: str, nsig: Union[None, float, List, Dict] = None, deltas_rel: Union[None, float] = None, allowNegativeSignals: bool = False)

Bases: object

Initialise.

Parameters

• dataObject – a smodels (combined)dataset or a list of theory predictions (for combination of analyses)

• nsig – signal yield, either as float or as list

• deltas_rel – relative error on signal. currently unused

AllowNegativeSignals

if True, negative values for the signal (mu) are allowed.

CLs(poi_test: float = 1.0, evaluationType: NllEvalType = NllEvalType.observed) → Optional[float]

compute CLs value for a given value of the poi

allowNegativeSignals

data

dataObject

dataType

deltas_sys

classmethod forAnalysesComb(theoryPredictions, deltas_rel)

get a statscomputer for combination of analyses :param theoryPredictions: list of TheoryPrediction objects :param deltas_rel: relative error for the signal :returns: a StatsComputer

classmethod forMultiBinSL(dataset, nsig, deltas_rel)

get a statscomputer for simplified likelihood combination.

Parameters

• dataset – CombinedDataSet object

• nsig – Number of signal events for each SR

Deltas_rel

Relative uncertainty for the signal

Returns

a StatsComputer

classmethod forPyhf(dataset, nsig, deltas_rel)

get a statscomputer for pyhf combination.

Parameters

• dataset – CombinedDataSet object

• nsig – Number of signal events for each SR

Deltas_rel

Relative uncertainty for the signal

Returns

a StatsComputer

classmethod forSingleBin(dataset, nsig, deltas_rel, lumi: Optional[Unum] = None)

get a statscomputer for an efficiency map (single bin).

Parameters

• dataset – DataSet object

• nsig – Number of signal events for each SR

Deltas_rel

Relative uncertainty for the signal

Returns

a StatsComputer

classmethod forTruncatedGaussian(theorypred, corr: float = 0.6)

get a statscomputer for truncated gaussians :param theorypred: TheoryPrediction object :param corr: correction factor: ULexp_mod = ULexp / (1. - corr*((ULobs-ULexp)/(ULobs+ULexp))) a factor of corr = 0.6 is proposed. :returns: a StatsComputer

getComputerAnalysesComb()

Create computer from a single bin :param nsig: signal yields.

getComputerMultiBinSL()

Create computer from a multi bin SL result

getComputerPyhf()

Create computer for a pyhf result

getComputerSingleBin(lumi)

Create computer from a single bin

getComputerTruncGaussian(**kwargs)

Create computer for truncated gaussians

get_five_values(evaluationType: NllEvalType, return_nll: bool = False, check_for_maxima: bool = False) → Dict

Return the Five Values: l(bsm), l(sm), muhat, l(muhat), sigma(mu_hat)

Parameters

check_for_maxima – if true, then check lmax against l(sm) and l(bsm) correct, if necessary

likelihood(poi_test: float, evaluationType: NllEvalType, return_nll: bool, asimov: Union[None, float] = None) → float

simple frontend to individual computers

likelihoodComputer

maximize_likelihood(evaluationType: NllEvalType, return_nll: bool = False) → dict

simple frontend to the individual computers, later spey :param return_nll: if True, return negative log likelihood :returns: Dictionary of llhd (llhd at mu_hat), muhat, sigma_mu (sigma of mu_hat), optionally also theta_hat

nsig

poi_upper_limit(evaluationType: NllEvalType, limit_on_xsec: bool = False, nSigma: int = 0) → Optional[Union[float, Unum]]

Simple frontend to the upperlimit computers, later to spey.poi_upper_limit

Parameters

limit_on_xsec – if True, then return the limit on the

cross section :param nSigma: the upper limit for central value (0), + 1 sigma, - 1 sigma, etc. For error bands.

restore(evaluationType)

SL only. transform the data to evaluationType or observed

transform(evaluationType)

SL only. transform the data to evaluationType or observed

upperLimitComputer

statistics.statsTools.getStatsComputerModule()

very single convenience function to centralize switching between our stats code and spey.

statistics.truncatedGaussians module

class statistics.truncatedGaussians.TruncatedGaussians(upperLimitOnMu: float, expectedUpperLimitOnMu: float, corr: Optional[float] = 0.6, cl=0.95)

Bases: object

likelihood computer based on the trunacated Gaussian approximation, see arXiv:1202.3415

Parameters

• upperLimitOnMu – observed upper limit on signal strength mu

• expectedUpperLimitOnMu – evaluationType upper limit on signal strength mu

• corr – correction factor: ULexp_mod = ULexp / (1. - corr*((ULobs-ULexp)/(ULobs+ULexp))) When comparing with likelihoods constructed from efficiency maps, a factor of corr = 0.6 has been found to result in the best approximations.

• cl – confidence level

cl

corr

denominator

expectedUpperLimitOnMu

likelihood(mu: Optional[float], return_nll: Optional[bool] = False, allowNegativeSignals: Optional[bool] = True, corr: Optional[float] = 0.6, evaluationType: NllEvalType = NllEvalType.observed) → Union[None, float]

return the likelihood, as a function of mu

Parameters

• mu – number of signal events, if None then mu = muhat

• return_nll – if True, return negative log likelihood

• allowNegativeSignals – if True, then allow muhat to become negative, else demand that muhat >= 0. In the presence of underfluctuations in the data, setting this to True results in more realistic approximate likelihoods.

Returns

likelihood (float)

lmax(return_nll: Optional[bool] = False, allowNegativeSignals: Optional[bool] = True, corr: Optional[float] = 0.6, evaluationType: NllEvalType = NllEvalType.observed) → Dict

Return the likelihood, as a function of mu

Parameters

• mu – number of signal events, if None then mu = muhat

• return_nll – if True, return negative log likelihood

• allowNegativeSignals – if True, then allow muhat to become negative, else demand that muhat >= 0. In the presence of underfluctuations in the data, setting this to True results in more realistic approximate likelihoods.

Returns

dictionary with likelihood (float), muhat, and sigma_mu

newCorrectionType = False

sigma_mu

upperLimitOnMu

Module contents