spey.backends.default_pdf.simple_pdf.MultivariateNormal

spey.backends.default_pdf.simple_pdf.MultivariateNormal#

class spey.backends.default_pdf.simple_pdf.MultivariateNormal(signal_yields: List[float] | Callable[[ndarray], ndarray], background_yields: List[float], data: List[int], covariance_matrix: List[List[float]] | callable, n_signal_parameters: int = 0, signal_parameter_bounds: List[Tuple[float | None, float | None]] | None = None)[source]#

Multivariate Normal likelihood with optional parameter-dependent covariance (default.multivariate_normal).

This backend models the likelihood as a single multivariate normal distribution over all bins simultaneously, capturing inter-bin correlations through an \(N \times N\) covariance matrix \(\Sigma\):

\[\mathcal{L}(\mu) = \frac{1}{\sqrt{(2\pi)^N \det\Sigma(\mu)}} \exp\!\left[-\frac{1}{2} \bigl(\boldsymbol{\lambda}(\mu) - \mathbf{n}^{\rm obs}\bigr)^\top \Sigma(\mu)^{-1} \bigl(\boldsymbol{\lambda}(\mu) - \mathbf{n}^{\rm obs}\bigr) \right],\]

where \(\lambda_i(\mu) = \mu n^s_i + n^b_i\) and \(\Sigma(\mu)\) may be constant or \(\mu\)-dependent (see below).

Callable covariance matrix

covariance_matrix can be a callable that takes the full parameter vector pars and returns an \(N \times N\) covariance matrix. This allows, for example, signal-strength-dependent uncertainties:

\[\Sigma(\mu) = \Sigma_{\rm bkg} + \mu^2\, \Sigma_{\rm sig}.\]

Added in version 0.1.9.

Changed in version 0.2.6: Callable covariance_matrix support added.

covariance_matrix can also take callable function as an input where function takes nuisance parameters as inputs and return a new covariance matrix as output.

Example — callable covariance matrix:

>>> import spey
>>> import numpy as np

>>> signal_yields = np.array([12.0, 15.0])
>>> background_yields = np.array([50.0, 48.0])
>>> data = np.array([36., 33.])
>>> covariance_matrix = np.array([[144.0, 13.0], [25.0, 256.0]])
>>> covariance_signal = np.array([[5.0, 1.0], [2.0, 3.0]])

>>> # ``pars`` include all the parameters of the model,
>>> # including the POI and any extra signal parameters
>>> def cov_matrix(pars: np.ndarray) -> np.ndarray:
>>>     return covariance_matrix + covariance_signal * pars[0]**2

>>> pdf_wrapper = spey.get_backend('default.multivariate_normal')
>>> model = pdf_wrapper(
...     signal_yields=signal_yields,
...     background_yields=background_yields,
...     data=data,
...     covariance_matrix=cov_matrix,
... )

Example — callable signal yields with additional parameters:

signal_yields can be a callable that accepts extra parameters beyond \(\mu\) and returns the per-bin signal yields as a np.ndarray. The number of extra parameters must be declared via n_signal_parameters so that the optimiser can construct the correct parameter vector. The first element of the parameter vector is always \(\mu\); the remaining elements (pars[1:]) are forwarded to the callable at each function evaluation.

>>> import spey
>>> import numpy as np

>>> background_yields = np.array([50.0, 48.0])
>>> data = np.array([36., 33.])
>>> covariance_matrix = np.array([[144.0, 13.0], [25.0, 256.0]])
>>> base_signal = np.array([12.0, 15.0])

>>> # signal_yields is a function of one extra parameter (signal_par_0)
>>> def signal_yields(extra_pars: np.ndarray) -> np.ndarray:
...     return base_signal * (1.0 + extra_pars[0])

>>> pdf_wrapper = spey.get_backend('default.multivariate_normal')
>>> model = pdf_wrapper(
...     signal_yields=signal_yields,
...     background_yields=background_yields,
...     data=data,
...     covariance_matrix=covariance_matrix,
...     n_signal_parameters=1,
... )

>>> muhat, nll = model.maximize_likelihood()

The resulting ModelConfig will contain two parameters: ["mu", "signal_par_0"].

Attention

The functional signal_yields is used to compute the expected signal yields in each bin. The optimiser only passes relevant parameter set to signal_yields function, excluding signal strength (i.e. the first parameter \(\mu\)). Callable covariance_matrix receives the full parameter vector, including the POI and any extra signal parameters.

Changed in version 0.2.6: The ability to input a callable covariance matrix has been added.

Changed in version 0.2.7: signal_yields now accepts a callable with signature (extra_pars: np.ndarray) -> np.ndarray in addition to a plain array. The companion argument n_signal_parameters (default 0) declares how many extra parameters the callable expects; they are appended to the parameter vector after \(\mu\) and their names are automatically set to signal_par_0, signal_par_1, …

Parameters:

  • signal_yields (List[float] | Callable[[np.ndarray], np.ndarray]) – signal yields per bin, or a callable that takes the extra signal parameters (pars[1:]) and returns the signal yields as a np.ndarray.

  • background_yields (List[float]) – background yields per bin.

  • data (List[int]) – observed data per bin.

  • covariance_matrix (List[List[float]] | callable) –

    covariance matrix (square matrix), or a callable that takes the full parameter vector and returns a covariance matrix.

  • n_signal_parameters (int, default 0) – number of additional free parameters to pass to a callable signal_yields. Has no effect when signal_yields is a plain array. When greater than zero the optimiser parameter vector is extended to [mu, signal_par_0, ..., signal_par_{n-1}].

  • signal_parameter_bounds (List[Tuple[Optional[float], Optional[float]]] | None) – Optimiser bounds for each extra signal parameter. Each entry is a (lower, upper) pair; use None for an unbounded side. When None, every extra signal parameter receives (None, None). Must have exactly n_signal_parameters entries when provided.

__init__(signal_yields: List[float] | Callable[[ndarray], ndarray], background_yields: List[float], data: List[int], covariance_matrix: List[List[float]] | callable, n_signal_parameters: int = 0, signal_parameter_bounds: List[Tuple[float | None, float | None]] | None = None)[source]#

Methods

__init__(signal_yields, background_yields, ...)

asimov_negative_loglikelihood([poi_test, ...])

Compute the profiled negative log-likelihood at fixed \(\mu\) on Asimov data.

combine(other, **kwargs)

Combine this statistical model with another backend instance.

config([allow_negative_signal, poi_upper_bound])

Model configuration.

expected_data(pars, **kwargs)

Compute the expected data vector \(\{\lambda_i(\mu)\}\) at the given parameter point.

get_hessian_logpdf_func([expected, data])

Return a callable that evaluates the Hessian of \(\ln\mathcal{L}(\mu)\).

get_logpdf_func([expected, data])

Return a callable that evaluates \(\ln\mathcal{L}(\mu)\).

get_objective_function([expected, data, do_grad])

Return the objective function \(-\ln\mathcal{L}(\mu)\) used by the optimiser.

get_sampler(pars)

Return a callable that draws pseudo-data from the main model.

minimize_asimov_negative_loglikelihood([...])

Find the global minimum of the negative log-likelihood on Asimov data (free fit).

minimize_negative_loglikelihood([expected, ...])

Find the global minimum of the negative log-likelihood (free fit).

negative_loglikelihood([poi_test, expected])

Compute the profiled negative log-likelihood at a fixed \(\mu\).

Attributes

covariance_matrix

n_signal_parameters

Number of signal parameters

author

Author of the backend

background_yields

data

is_alive

Returns True if at least one bin has non-zero signal yield.

main_model

Main model distribution — Poisson (or Gaussian) term of the likelihood.

name

Name of the backend

signal_yields

spey_requires

Spey version required for the backend

version

Version of the backend

constraints

Constraints to be used during optimisation process