derivkit package#

Subpackages#

Submodules#

Module contents#

Provides all derivkit methods.

class derivkit.CalculusKit(function: Callable[[Sequence[float] | np.ndarray], float | NDArray[np.floating]], x0: Sequence[float] | np.ndarray, thread_safe: bool = False, thread_lock: threading.RLock | None = None)#

Bases: object

Provides access to gradient, Jacobian, Hessian, and higher-order derivative tensors.

Initialises class with function and expansion point.

Parameters:

function – The function to be differentiated. Accepts a 1D array-like. Must return either a scalar or a 1D array, depending on the requested derivative tensor.
x0 – Point at which to evaluate derivatives (shape (p,)) for p input parameters.
thread_safe – If True, serialize calls to function using a lock. This prevents concurrent evaluation within a single Python process (i.e. thread-based parallelism). It does not synchronize across multiple processes.
thread_lock – Optional lock object to use when thread_safe=True. If None, an internal lock is created. Ignored if thread_safe=False.

gradient(*args, **kwargs) → NDArray[floating]#

Returns the gradient of a scalar-valued function.

This is a wrapper around derivkit.calculus.build_gradient(), with the function and theta0 arguments fixed to the values provided at initialization of CalculusKit. Any additional positional or keyword arguments are forwarded to derivkit.calculus.build_gradient().

Refer to the documentation of derivkit.calculus.build_gradient() for available options.

hessian(*args, **kwargs) → NDArray[floating]#

Returns the Hessian of a scalar-valued function.

This is a wrapper around derivkit.calculus.build_hessian(), with the function and theta0 arguments fixed to the values provided at initialization of CalculusKit. Any additional positional or keyword arguments are forwarded to derivkit.calculus.build_hessian().

Refer to the documentation of derivkit.calculus.build_hessian() for available options.

hessian_diag(*args, **kwargs) → NDArray[floating]#

Returns the diagonal of the Hessian of a scalar-valued function.

This is a wrapper around derivkit.calculus.build_hessian_diag(), with the function and theta0 arguments fixed to the values provided at initialization of CalculusKit. Any additional positional or keyword arguments are forwarded to derivkit.calculus.build_hessian_diag().

Refer to the documentation of derivkit.calculus.build_hessian_diag() for available options.

hyper_hessian(*args, **kwargs) → NDArray[floating]#

Returns a higher-order derivative tensor of a function.

This is a wrapper around derivkit.calculus.build_hyper_hessian(), with the function and theta0 arguments fixed to the values provided at initialization of CalculusKit. Any additional positional or keyword arguments are forwarded to derivkit.calculus.build_hyper_hessian().

Refer to the documentation of derivkit.calculus.build_hyper_hessian() for available options.

jacobian(*args, **kwargs) → NDArray[floating]#

Returns the Jacobian of a vector-valued function.

This is a wrapper around derivkit.calculus.build_jacobian(), with the function and theta0 arguments fixed to the values provided at initialization of CalculusKit. Any additional positional or keyword arguments are forwarded to derivkit.calculus.build_jacobian().

Refer to the documentation of derivkit.calculus.build_jacobian() for available options.

Bases: object

Unified interface for computing numerical derivatives.

The class provides a simple way to evaluate derivatives using any of DerivKit’s available backends (e.g., adaptive fit or finite difference). By default, the local polynomial fitting method is used.

You can supply either a function and x0, or tabulated tab_x/tab_y and x0 in case you want to differentiate a tabulated function. The chosen backend is invoked when you call the .differentiate() method.

The input function can be wrapped with an input-based LRU cache for numerical derivative backends. Automatic differentiation backends use the raw function to preserve traceability.

Example

>>> import numpy as np
>>> from derivkit import DerivativeKit
>>> dk = DerivativeKit(np.cos, x0=1.0)
>>> deriv = dk.differentiate(order=1)  # uses the default method

function#: The callable used for derivative evaluation.

x0#: The point or points at which the derivative is evaluated.

default_method#: The numerical differentiation method used when no method is specified. Initially "local_polynomial".

Initializes the DerivativeKit with a target function and expansion point.

Parameters:

function – The function to be differentiated. Must accept a single float and return a scalar or array-like output.
x0 – Point or array of points at which to evaluate the derivative.
tab_x – Optional tabulated x values for creating a tabulated_model.one_d.Tabulated1DModel.
tab_y – Optional tabulated y values for creating a tabulated_model.one_d.Tabulated1DModel.
use_input_cache – If True, wrap the input function with an input-based LRU cache. Exact input values are cached by default.
cache_number_decimal_places – Decimal places used when constructing cache keys. If None, exact input values are used.
cache_maxsize – Maximum size of the input cache.
cache_copy – If True, return copies of cached NumPy array outputs.

cache_clear() → None#: Clears the cache.

cache_info() → Any#: Returns cache info about the cache.

differentiate(*, method: str | None = None, **kwargs: Any) → Any#

Compute derivatives using the chosen method.

Forwards all keyword arguments to the engine’s .differentiate().

Parameters:

method – Numerical differentiation method or alias (e.g., "adaptive", "finite"). If None, uses self.default_method.
**kwargs – Passed through to the chosen engine.

Returns:

The derivative result from the underlying engine.

If x0 is a single value, returns the usual derivative output.

If x0 is an array of points, returns an array where the first dimension indexes the points in x0. For example, if you pass 5 points and each derivative has shape (2, 3), the result has shape (5, 2, 3).

Notes

Thread-level parallelism across derivative evaluations can be controlled by passing n_workers via **kwargs. Note that this does not launch separate Python processes. All work occurs within a single process using worker threads.

Raises:: ValueError – If method is not recognized.

class derivkit.ForecastKit(function: Callable[[Sequence[float] | np.ndarray], np.ndarray] | None, theta0: Sequence[float] | np.ndarray, cov: np.ndarray | Callable[[np.ndarray], np.ndarray], *, use_input_cache: bool = True, cache_number_decimal_places: int | None = None, cache_maxsize: int | None = 4096, cache_copy: bool = True, thread_safe: bool = False, thread_lock: threading.RLock | None = None)#

Bases: object

Provides access to forecasting workflows.

Initialises the ForecastKit with model, fiducials, and covariance.

As an optimization the class can cache function values, which speeds up repeated function calls and their associated derivatives. In this case the function will truncate the number of decimals of the input parameters.

Parameters:

function – Callable returning the model mean vector \(\mu(\theta)\). May be None if you only plan to use covariance-only workflows (e.g. generalized Fisher with term="cov"). Required for fisher(), dali(), and fisher_bias().
theta0 – Fiducial parameter values of shape (p, ) where p is the number of parameters.
cov –
Covariance specification. Supported forms are:
- cov=C0: fixed covariance matrix \(C(\theta_0)\) with shape (n_obs, n_obs), where n_obs is the number of observables.
- cov=cov_fn: callable with cov_fn(theta) returning the covariance matrix \(C(\theta)\) evaluated at the parameter vector theta, with shape (n_obs, n_obs). The covariance at theta0 is evaluated once and cached.
use_input_cache – If True, wrap the input function with an input-based LRU cache. Exact input values are cached by default.
cache_number_decimal_places – Decimal places used when constructing cache keys. If None, exact input values are used.
cache_maxsize – Maximum size of the input cache.
cache_copy – If True, return copies of cached NumPy array outputs.
thread_safe – If True, serialize calls to function using a lock. This prevents concurrent evaluation within a single Python process (i.e. thread-based parallelism). It does not synchronize across multiple processes.
thread_lock – Optional lock object to use when thread_safe=True. If None, an internal lock is created.

dali(*, method: str | None = None, forecast_order: int = 2, symmetrize: bool = True, n_workers: int = 1, return_derivatives: bool = False, **dk_kwargs: Any) → dict[int, tuple[NDArray[float64], ...]] | tuple[dict[int, tuple[NDArray[float64], ...]], dict[str, Any]]#

Builds the DALI expansion for the given model up to the given order.

Parameters:

method – Method name or alias (e.g., "adaptive", "finite"). If None, the derivkit.derivative_kit.DerivativeKit default is used.
forecast_order – The requested order of the forecast. Currently supported values and their meaning are given in derivkit.forecasting.forecast_core.SUPPORTED_FORECAST_ORDERS.
symmetrize – Flag to force symmetrization across all DALI tensor axes.
n_workers – Number of workers for per-parameter parallelization/threads. Default 1 (serial). Inner batch evaluation is kept serial to avoid oversubscription.
return_derivatives – When set to True, returns the underlying derivatives.
**dk_kwargs – Additional keyword arguments passed to derivkit.calculus_kit.CalculusKit.

Returns:

A dict mapping order -> multiplet for all order = 1..forecast_order.

For each forecast order k, the returned multiplet contains the tensors introduced at that order. Concretely:

order 1: (F_{(1,1)},) (Fisher matrix)
order 2: (D_{(2,1)}, D_{(2,2)})
order 3: (T_{(3,1)}, T_{(3,2)}, T_{(3,3)})
order 4: (Qa_{(4,1)}, Qa_{(4,2)}, Qa_{(4,3)}, Qa_{(4,4)})

Here the tensors at order k are contractions of the k-th order model derivative with derivatives of orders 1 through k via the inverse covariance.

Each tensor axis has length p = len(self.theta0). The additional tensors at order k have parameter-axis ranks from k+1 through 2*k.

If return_derivatives is set to True: a tuple with the above dict as the first element, and the dict of derivatives as the second element.

dali_bias_tensor(*, delta_nu: Sequence[float] | NDArray[floating] | Sequence[Sequence[float]], bias_order: int = 3, method: str | None = None, n_workers: int = 1, derivatives: dict[int, NDArray[float64]] | None = None, **dk_kwargs: Any) → dict[int, NDArray[float64]]#

Builds the systematic mismatch tensors entering the DALI bias expansion.

The tensors describe how the data-model mismatch delta_nu couples to successive derivatives of the model.

Parameters:

delta_nu – Difference between two data vectors, which may represent a systematic mismatch or a difference between model predictions. Accepts a 1D array or a column vector of shape (n_observables, 1). The input is flattened internally.
bias_order – Highest mismatch tensor order to compute. Supported values are given in derivkit.forecasting.forecast_core.SUPPORTED_DERIVATIVE_ORDERS.
method – Numerical differentiation method or alias (e.g., "adaptive", "finite"). If None, uses derivkit.derivative_kit.DerivativeKit.default_method.
n_workers – Number of workers for per-parameter parallelization/threads. Default 1 (serial).
derivatives – Optional dictionary of model derivatives keyed by order. Existing derivatives are reused. Missing orders are computed and added to the supplied dictionary in place, modifying the input. If None, an internal dictionary is used.
**dk_kwargs – Additional keyword arguments passed to derivkit.derivative_kit.DerivativeKit.differentiate().

Returns:

A dictionary mapping derivative orders 1 up to and including bias_order to their corresponding mismatch tensors. Zeroth-order terms are omitted because they are independent of the parameter displacement.

Raises:

TypeError – If bias_order is not an integer.
ValueError – If bias_order is unsupported, theta0 is empty, or the data mismatch does not match the number of observables.

dali_map_shift(*, delta_nu: Sequence[float] | NDArray[floating] | Sequence[Sequence[float]], expansion_order: int = 3, return_metadata: bool = False, method: str | None = None, n_workers: int = 1, derivatives: dict[int, NDArray[float64]] | None = None, **dk_kwargs: Any)#

Predict how the posterior maximum moves when a systematic is introduced.

Parameters:

delta_nu – Actual difference between the systematically affected data vector and the fiducial model prediction. Its size is already included; no separate amplitude parameter is needed.
expansion_order – Number of perturbative MAP response orders (1 to 3).
return_metadata – If set to True, return the shift and a dictionary of explanatory results; otherwise return only the shift.
method – Numerical differentiation method or alias (e.g., "adaptive", "finite"). If None, uses derivkit.derivative_kit.DerivativeKit.default_method.
n_workers – Number of workers for per-parameter parallelization/threads. Default 1 (serial).
derivatives – Optional dictionary of model derivatives keyed by order. Existing derivatives are reused. Missing orders are computed and added to the supplied dictionary in place, modifying the input. If None, an internal dictionary is used.
**dk_kwargs – Additional keyword arguments passed to derivkit.derivative_kit.DerivativeKit.differentiate().

Returns:

If return_metadata=False, an array giving the MAP systematic shift: MAP with the mismatch minus MAP without the mismatch. Since the baseline MAP is assumed to be theta0, this is also the MAP displacement from the fiducial parameters.

If return_metadata=True, returns (shift, metadata). The metadata dictionary contains:

estimator: Predicted MAP parameter values with the mismatch.: This is theta0 + shift.

baseline: MAP without the mismatch, assumed to be theta0.

displacement: MAP with the mismatch minus theta0.: This equals shift because the baseline MAP is assumed to be theta0.
terms: Dictionary keyed by response order. Entry 1 is the: leading Fisher bias; entries 2 and 3 are the successive corrections. Their sum is the returned shift.
gaussian_terms: Terms obtained if the unbiased posterior is: approximated by its quadratic (Gaussian) Fisher form, while retaining the mismatch potential up to and including the requested order.
non_gaussian_terms: Differences between the DALI and Gaussian: terms at each order. They isolate corrections caused by the non-quadratic reference posterior within this expansion.

gaussian_contribution: Sum of gaussian_terms.

non_gaussian_contribution: Sum of non_gaussian_terms.: Adding the two contributions reproduces the MAP shift.

These Gaussian/non-Gaussian labels refer to the reference posterior, not to whether a term is linear or nonlinear in systematic amplitude.

dali_mean_response(*, delta_nu: Sequence[float] | NDArray[floating] | Sequence[Sequence[float]], theta: NDArray[float64], weights: NDArray[float64] | None = None, forecast_order: int = 3, response_order: int = 4, return_metadata: bool = False, method: str | None = None, n_workers: int = 1, derivatives: dict[int, NDArray[float64]] | None = None, **dk_kwargs: Any)#

Predict the posterior mean shift from its systematic response.

Expands the mean response around an unbiased reference posterior represented by samples or integration nodes.

Parameters:

delta_nu – Full systematic difference in the data vector.
theta – Samples or integration nodes representing the unbiased reference posterior, with shape (n_samples, n_parameters).
weights – Weights representing the reference posterior at the supplied points. If None, equal weights are used, appropriate for equally weighted posterior samples. For integration nodes, the weights must include the reference posterior density and integration measure. Weights are normalized internally.
forecast_order – Highest model derivative order included in the systematic mismatch potential. Supported values are given in derivkit.forecasting.forecast_core.SUPPORTED_FORECAST_ORDERS.
response_order – Highest order retained in the systematic amplitude expansion of the posterior mean. Must be 1, 2, 3, or 4.
return_metadata – If True, return the shift and explanatory metadata.
method – Numerical differentiation method or alias (e.g., "adaptive", "finite"). If None, uses derivkit.derivative_kit.DerivativeKit.default_method
n_workers – Number of workers for per-parameter parallelization/threads. Default 1 (serial).
derivatives – Optional dictionary of model derivatives keyed by order. Existing derivatives are reused. Missing orders are computed and added to the supplied dictionary in place, modifying the input. If None, an internal dictionary is used.
**dk_kwargs – Additional keyword arguments passed to derivkit.derivative_kit.DerivativeKit.differentiate().

Returns:

An array containing the predicted systematic shift of the posterior mean, calculated by summing the response terms through response_order. The shift is the posterior mean with the mismatch minus the posterior mean without it.

If return_metadata=True, returns (shift, metadata) instead. The metadata dictionary contains:

estimator: Predicted posterior mean with the mismatch.: This is the reference posterior mean plus the predicted shift.
baseline: Posterior mean without the mismatch, calculated from: the supplied samples.
displacement: Predicted posterior mean with the mismatch minus: the fiducial parameters theta0.
terms: Dictionary mapping each response order to its contribution: to the predicted posterior mean shift for the supplied mismatch. Their sum equals the returned shift.

Raises:

TypeError – If forecast_order or response_order is not an integer.
ValueError – If either order is unsupported, or if the supplied reference samples have invalid shapes.

dali_mean_shift(*, delta_nu: Sequence[float] | NDArray[floating] | Sequence[Sequence[float]], return_metadata: bool = False, method: str | None = None, n_workers: int = 1, derivatives: dict[int, NDArray[float64]] | None = None, **dk_kwargs: Any)#

Predict the systematic shift of the posterior mean analytically.

Uses third-order DALI tensors to calculate the perturbative MAP shift and the change in the posterior mean relative to the MAP.

Parameters:

delta_nu – Difference between the systematically affected data vector and the fiducial model prediction, including its magnitude.
return_metadata – If True, return the shift and explanatory metadata.
method – Numerical differentiation method or alias (e.g., "adaptive", "finite"). If None, uses derivkit.derivative_kit.DerivativeKit.default_method.
n_workers – Number of workers for per-parameter parallelization/threads. Default 1 (serial).
derivatives – Optional dictionary of model derivatives keyed by order. Existing derivatives are reused. Missing orders are computed and added to the supplied dictionary in place, modifying the input. If None, an internal dictionary is used.
**dk_kwargs – Additional keyword arguments passed to derivkit.derivative_kit.DerivativeKit.differentiate().

Returns:

An array containing the predicted systematic shift of the posterior mean, defined as the posterior mean with the mismatch minus the posterior mean without it.

If return_metadata=True, returns (shift, metadata) instead. metadata is a dictionary containing:

estimator: Predicted posterior mean with the mismatch.

baseline: Predicted posterior mean without the mismatch.: This need not equal theta0.
displacement: Predicted posterior mean with the mismatch minus: the fiducial parameters theta0.

map_shift: Predicted systematic shift of the posterior MAP.

delta_chi2_dali(*, theta: ndarray, dali: dict[int, tuple[ndarray, ...]], forecast_order: int | None = 2) → float#

Computes a displacement chi-squared under the DALI approximation.

This evaluates a scalar delta_chi2 from the displacement d = theta - theta0 using the provided DALI tensors.

The expansion point is taken from ForecastKit.theta0.

Parameters:

theta – Evaluation point in parameter space with shape (p,).
dali – DALI tensors as returned by ForecastKit.dali().
forecast_order – Maximum order of the DALI expansion to include. If None, the highest order available in dali is used.

Returns:

Scalar delta chi-squared value.

delta_chi2_fisher(*, theta: ndarray, fisher: ndarray) → float#

Computes a displacement chi-squared under the Fisher approximation.

This evaluates the standard quadratic form

delta_chi2 = (theta - theta0)^T @ F @ (theta - theta0)

using the provided Fisher matrix and the stored expansion point ForecastKit.theta0.

Parameters:

theta – Evaluation point in parameter space with shape (p,).
fisher – Fisher matrix with shape (p, p).

Returns:

Scalar delta chi-squared value.

delta_nu(data_unbiased: ndarray, data_biased: ndarray) → ndarray#

Computes the difference between two data vectors.

This helper is used in Fisher-bias calculations and any other workflow where two data vectors are compared: it takes a pair of vectors (for example, a version with a systematic and one without) and returns their difference as a 1D array whose length matches the number of observables implied by cov. It works with both 1D inputs and 2D arrays (for example, correlation-by-ell) and flattens 2D inputs using NumPy’s row-major (“C”) order, which is the standard convention throughout the DerivKit package.

Parameters:

data_unbiased – Reference data vector without the systematic. Can be 1D or 2D. If 1D, it must follow the NumPy’s row-major (“C”) flattening convention used throughout the package.
data_biased – Data vector that includes the systematic effect. Can be 1D or 2D. If 1D, it must follow the NumPy’s row-major (“C”) flattening convention used throughout the package.

Returns:

A 1D NumPy array of length n_observables representing the mismatch between the two input data vectors. This is simply the element-wise difference between the input with systematic and the input without systematic, flattened if necessary to match the expected observable ordering.

fisher(*, method: str | None = None, n_workers: int = 1, return_derivatives: bool = False, **dk_kwargs: Any) → ndarray | tuple[ndarray, ndarray]#

Computes the Fisher information matrix for a given model and covariance.

Parameters:

method – Derivative method name or alias (e.g., "adaptive", "finite"). If None, the derivkit.derivative_kit.DerivativeKit default is used.
n_workers – Number of workers for per-parameter parallelisation. Default is 1 (serial).
return_derivatives – When set to True, returns the underlying derivatives.
**dk_kwargs – Additional keyword arguments forwarded to derivkit.derivative_kit.DerivativeKit.differentiate().

Returns:

Returns the Fisher matrix with shape (n_parameters, n_parameters). If return_derivatives is set to True, returns a tuple with the Fisher matrix as its first element, while the second element is the the underlying Jacobian.

fisher_bias(*, fisher_matrix: ndarray, delta_nu: ndarray, method: str | None = None, n_workers: int = 1, rcond: float = 1e-12, **dk_kwargs: Any) → tuple[ndarray, ndarray]#

Estimates parameter bias using the stored model, expansion point, and covariance.

This function takes a model, an expansion point, a covariance matrix, a Fisher matrix, and a data-vector difference delta_nu and maps that difference into parameter space. A common use case is the classic “Fisher bias” setup, where one asks how a systematic-induced change in the data would shift inferred parameters.

Internally, the function evaluates the model response at the expansion point and uses the covariance and Fisher matrix to compute both the parameter-space bias vector and the corresponding shifts. See https://arxiv.org/abs/0710.5171 for details.

Parameters:

fisher_matrix – Square matrix describing information about the parameters. Its shape must be (p, p), where p is the number of parameters.
delta_nu – Difference between a biased and an unbiased data vector, for example \(\Delta\nu = \nu_{\mathrm{biased}} - \nu_{\mathrm{unbiased}}\). Accepts a 1D array of length n or a 2D array that will be flattened in row-major order (“C”) to length n, where n is the number of observables. If supplied as a 1D array, it must already follow the same row-major (“C”) flattening convention used throughout the package.
n_workers – Number of workers used by the internal derivative routine when forming the Jacobian.
method – Method name or alias (e.g., "adaptive", "finite"). If None, the derivkit.derivative_kit.DerivativeKit default is used.
rcond – Regularization cutoff for pseudoinverse. Default is 1e-12.
**dk_kwargs – Additional keyword arguments passed to derivkit.derivative_kit.DerivativeKit.differentiate().

Returns:

A tuple (bias_vec, delta_theta) of 1D arrays with length p, where bias_vec is the parameter-space bias vector and delta_theta are the corresponding parameter shifts.

gaussian_fisher(*, method: str | None = None, n_workers: int = 1, rcond: float = 1e-12, symmetrize_dcov: bool = True, return_derivatives: bool = False, **dk_kwargs: Any) → ndarray | tuple[ndarray, dict[str, ndarray | None]]#

Computes the generalized Fisher matrix for parameter-dependent mean and/or covariance.

This function computes the generalized Fisher matrix for a Gaussian likelihood with parameter-dependent mean and/or covariance. Uses derivkit.forecasting.fisher_gaussian.build_gaussian_fisher_matrix().

Parameters:

method – Derivative method name or alias (e.g., "adaptive", "finite").
n_workers – Number of workers for per-parameter parallelisation.
rcond – Regularization cutoff for pseudoinverse fallback in linear solves.
symmetrize_dcov – If True, symmetrize each covariance derivative via \(\tfrac{1}{2}(C_{,i} + C_{,i}^{\mathsf{T}})\).
return_derivatives – If True, return the underlying Jacobian.
**dk_kwargs – Forwarded to the internal derivative calls.

Returns:

Fisher matrix with shape (p, p) with p the number of parameters.

If return_metadata is True, return a tuple containing the Fisher matrix as the first element and the Jacobian as the second element.

getdist_dali_emcee(*, dali: dict[int, tuple[ndarray, ...]], bias: dict[int, ndarray[floating]] | None = None, names: Sequence[str], labels: Sequence[str], **kwargs: Any)#

Returns GetDist getdist.MCSamples from emcee sampling of a DALI posterior.

This is a thin wrapper around derivkit.forecasting.getdist_dali_samples.dali_to_getdist_emcee() that fixes the expansion point to self.theta0.

Parameters:

dali – DALI tensors as returned by ForecastKit.dali().
bias – Optional dictionary of systematic bias tensors indexed by expansion order, as returned by derivkit.forecasting.build_dali_bias_tensor(). These contribute to the target log posterior.
names – Parameter names for GetDist (length p).
labels – Parameter labels for GetDist (length p).
**kwargs – Forwarded to derivkit.forecasting.getdist_dali_samples.dali_to_getdist_emcee() (e.g. forecast_order, n_steps, burn, thin, n_walkers, init_scale, seed, prior_terms, prior_bounds, logprior, sampler_bounds, label).

Returns:

A getdist.MCSamples containing MCMC chains.

getdist_dali_importance(*, dali: dict[int, tuple[ndarray, ...]], bias: dict[int, ndarray[floating]] | None = None, names: Sequence[str], labels: Sequence[str], **kwargs: Any)#

Returns GetDist getdist.MCSamples for a DALI posterior via importance sampling.

This is a thin wrapper around derivkit.forecasting.getdist_dali_samples.dali_to_getdist_importance() that fixes the expansion point to self.theta0.

Parameters:

dali – DALI tensors as returned by ForecastKit.dali().
bias – Optional dictionary of systematic bias tensors indexed by expansion order, as returned by derivkit.forecasting.build_dali_bias_tensor(). These contribute to the target log posterior.
names – Parameter names for GetDist (length p).
labels – Parameter labels for GetDist (length p).
**kwargs – Forwarded to derivkit.forecasting.getdist_dali_samples.dali_to_getdist_importance() (e.g. forecast_order, n_samples, kernel_scale, seed, prior_terms, prior_bounds, logprior, sampler_bounds, label).

Returns:

A getdist.MCSamples with importance weights.

getdist_fisher_gaussian(*, fisher: ndarray, bias: ndarray | None = None, names: Sequence[str] | None = None, labels: Sequence[str] | None = None, **kwargs: Any)#

Converts a Fisher Gaussian into a GetDist getdist.gaussian_mixtures.GaussianND.

This is a thin wrapper around derivkit.forecasting.getdist_fisher_samples.fisher_to_getdist_gaussiannd() that fixes the mean to the stored expansion point self.theta0.

Parameters:

fisher – Fisher matrix with shape (p, p) evaluated at self.theta0.
bias – Optional Fisher systematic bias vector of shape (p,). Shifts the Gaussian mean by covariance @ bias.
names – Optional parameter names (length p).
labels – Optional parameter labels (length p).
**kwargs – Forwarded to derivkit.forecasting.getdist_fisher_samples.fisher_to_getdist_gaussiannd() (e.g. label, rcond).

Returns:

A getdist.gaussian_mixtures.GaussianND with mean self.theta0 and covariance given by the (pseudo-)inverse Fisher matrix.

getdist_fisher_samples(*, fisher: ndarray, bias: ndarray | None = None, names: Sequence[str], labels: Sequence[str], **kwargs: Any)#

Draws GetDist getdist.MCSamples from the Fisher Gaussian at self.theta0.

This is a thin wrapper around derivkit.forecasting.getdist_fisher_samples.fisher_to_getdist_samples() that fixes the sampling center to the stored expansion point self.theta0.

Parameters:

fisher – Fisher matrix with shape (p, p) evaluated at self.theta0.
bias – Optional Fisher systematic bias vector of shape (p,). Shifts the Gaussian sampling mean by pinv(fisher) @ bias and contributes to the stored negative log posterior.
names – Parameter names for GetDist (length p).
labels – Parameter labels for GetDist (length p).
**kwargs – Forwarded to derivkit.forecasting.getdist_fisher_samples.fisher_to_getdist_samples() (e.g. n_samples, seed, kernel_scale, prior_terms, prior_bounds, logprior, hard_bounds, store_loglikes, label).

Returns:

A getdist.MCSamples object containing samples drawn from the Fisher Gaussian.

laplace_approximation(*, neg_logposterior: Callable[[ndarray], float], theta_map: Sequence[float] | ndarray | None = None, method: str | None = None, n_workers: int = 1, ensure_spd: bool = True, rcond: float = 1e-12, **dk_kwargs: Any) → dict[str, Any]#

Computes a Laplace (Gaussian) approximation around theta_map.

The Laplace approximation replaces the posterior near its peak with a Gaussian. It does this by measuring the local curvature of the negative log-posterior using its Hessian at theta_map. The Hessian acts like a local precision matrix, and its inverse is the approximate covariance.

If theta_map is not provided, this uses the stored expansion point self.theta0.

Parameters:

neg_logposterior – Callable that accepts a 1D float64 parameter vector and returns a scalar negative log-posterior value.
theta_map – Expansion point for the approximation. This is often the maximum a posteriori estimate (MAP). If None, uses self.theta0.
method – Derivative method name/alias forwarded to the Hessian builder.
n_workers – Outer parallelism forwarded to Hessian construction.
ensure_spd – If True, attempt to regularize the Hessian to be symmetric positive definite (SPD) by adding diagonal jitter.
rcond – Cutoff for small singular values used by the pseudoinverse fallback when computing the covariance.
**dk_kwargs – Additional keyword arguments forwarded to derivkit.derivative_kit.DerivativeKit.differentiate().

Returns:

Dictionary with the Laplace approximation outputs (theta_map, neg_logposterior_at_map, hessian, cov, and jitter).

laplace_covariance(hessian: ndarray, *, rcond: float = 1e-12) → ndarray#

Computes the Laplace covariance matrix from a Hessian.

In the Laplace (Gaussian) approximation, the Hessian of the negative log-posterior at the expansion point acts like a local precision matrix. The approximate posterior covariance is the matrix inverse of that Hessian.

Parameters:

hessian – 2D square Hessian matrix.
rcond – Cutoff for small singular values used by the pseudoinverse fallback.

Returns:

A 2D symmetric covariance matrix with the same shape as hessian.

laplace_hessian(*, neg_logposterior: Callable[[ndarray], float], theta_map: Sequence[float] | ndarray | None = None, method: str | None = None, n_workers: int = 1, **dk_kwargs: Any) → ndarray#

Computes the Hessian of the negative log-posterior at theta_map.

The Hessian at theta_map measures the local curvature of the posterior peak. In the Laplace approximation, this Hessian plays the role of a local precision matrix, and its inverse provides a fast Gaussian estimate of parameter uncertainties and correlations.

If theta_map is not provided, this uses the stored expansion point self.theta0.

Parameters:

neg_logposterior – Callable returning the scalar negative log-posterior.
theta_map – Point where the curvature is evaluated (typically the MAP). If None, uses self.theta0.
method – Derivative method name/alias forwarded to the calculus machinery.
n_workers – Outer parallelism forwarded to Hessian construction.
**dk_kwargs – Additional keyword arguments forwarded to derivkit.derivative_kit.DerivativeKit.differentiate().

Returns:

A symmetric 2D array with shape (p, p) giving the Hessian of neg_logposterior evaluated at theta_map.

Computes the log posterior (up to a constant) under the DALI approximation.

If no prior is provided, this returns the DALI log-likelihood expansion with a flat prior and no hard cutoffs. Priors may be provided either as a pre-built logprior(theta) callable or as a lightweight prior specification via prior_terms and/or prior_bounds.

The expansion point is taken from the stored self.theta0.

Parameters:

theta – Evaluation point in parameter space with shape (p,).
dali – DALI tensors as returned by ForecastKit.dali().
bias – Optional systematic-bias tensors mapping expansion order to parameter-space tensors. Adds the factorial-weighted mismatch expansion through forecast_order. If None, evaluates the unbiased DALI posterior.
forecast_order – Maximum order of the DALI expansion to include. If None, the highest order available in dali is used.
prior_terms – Prior term specification passed to the underlying prior builder. Use this only if logprior is not provided.
prior_bounds – Global hard bounds passed to the underlying prior builder. Use this only if logprior is not provided.
logprior – Optional custom log-prior callable. If it returns a non-finite value, the posterior is treated as zero at that point and the function returns -np.inf. Cannot be used together with prior_terms or prior_bounds.

Returns:

Scalar log posterior value, defined up to an additive constant.

Computes the log posterior under the Fisher approximation.

The returned value is defined up to an additive constant in log space. This corresponds to an overall multiplicative normalization of the posterior density in probability space.

If no prior is provided, this returns the Fisher log-likelihood expansion with a flat prior and no hard cutoffs. Priors may be provided either as a pre-built logprior(theta) callable or as a lightweight prior specification via prior_terms and/or prior_bounds.

The expansion point is taken from the stored self.theta0.

Parameters:

theta – Evaluation point in parameter space with shape (p,).
fisher – Fisher matrix with shape (p, p).
bias – Optional linear systematic-bias vector of shape (p,). If provided, adds bias @ (theta - theta0) to the Fisher log posterior.
prior_terms – Prior term specification passed to the underlying prior builder. Use this only if logprior is not provided.
prior_bounds – Global hard bounds passed to the underlying prior builder. Use this only if logprior is not provided.
logprior – Optional custom log-prior callable. If it returns a non-finite value, the posterior is treated as zero at that point and the function returns -np.inf. Cannot be used together with prior_terms or prior_bounds.

Returns:

Scalar log posterior value, defined up to an additive constant.

negative_logposterior(theta: Sequence[float] | ndarray, *, logposterior: Callable[[ndarray], float]) → float#

Computes the negative log-posterior at theta.

This converts a log-posterior callable into the objective used by MAP estimation and curvature-based methods. It simply returns -logposterior(theta) and validates that the result is finite.

Parameters:

theta – 1D array-like parameter vector.
logposterior – Callable that accepts a 1D float64 array and returns a scalar float.

Returns:

Negative log-posterior value as a float.

xy_fisher(*, x0: ndarray, mu_xy: Callable[[ndarray, ndarray], ndarray], cov_xy: ndarray, method: str | None = None, n_workers: int = 1, rcond: float = 1e-12, symmetrize_dcov: bool = True, check_cyy_consistency: bool = True, atol: float = 0.0, rtol: float = 1e-12, return_derivatives: bool = False, **dk_kwargs: Any) → ndarray | tuple[ndarray, dict[str, ndarray]]#

Computes the X–Y Gaussian Fisher matrix (noisy inputs and outputs).

This implements the X–Y Gaussian Fisher formalism where the measured inputs and outputs are both noisy and may be correlated. Input uncertainty is propagated into an effective output covariance via a local linearization of mu_xy(x, theta) with respect to x at the measured inputs x0.

Parameters:

x0 – Measured inputs with shape (n_x,).
mu_xy – Model mean callable mu_xy(x, theta) -> y.
cov_xy – Full covariance for the stacked vector [x, y] with shape (n_x + n_y, n_x + n_y). Note that this is the joint covariance of the input and output measurements, ordered as [x0, y0] (inputs first, then outputs), so it contains the blocks Cxx, Cxy, and Cyy.
method – Derivative method name/alias forwarded to the derivative backend.
n_workers – Number of workers for derivative evaluations.
rcond – Regularization cutoff used in linear solves / pseudoinverse fallback.
symmetrize_dcov – If True, symmetrize each covariance derivative.
check_cyy_consistency – If True and this ForecastKit was initialized with a fixed output covariance cov=self.cov0 (i.e. \(C_{yy}\)), verify that it matches the \(C_{yy}\) block implied by cov_xy.
atol – Absolute tolerance for the \(C_{yy}\) consistency check.
rtol – Relative tolerance for the \(C_{yy}\) consistency check.
return_derivatives – When set to True, returns the derivatives used to build the Fisher matrix.
**dk_kwargs – Forwarded to the derivative backend.

Returns:

Fisher matrix with shape (p, p), where p is the number of parameters.

When return_derivatives is set to True, returns a tuple containing the Fisher matrix as the first element, while the derivatives are the second element.

Raises:

ValueError – If check_cyy_consistency=True and self.cov0 is compatible with a \(C_{yy}\) block but differs from the \(C_{yy}\) block extracted from cov_xy beyond rtol/atol.

class derivkit.LikelihoodKit(data: ArrayLike, model_parameters: ArrayLike)#

Bases: object

High-level interface for Gaussian and Poissonian likelihoods.

The class stores data and model_parameters and provides methods to evaluate the corresponding likelihoods.

Initialises the likelihoods object.

Parameters:

data – Observed data values. The expected shape depends on the particular likelihoods. For the Gaussian likelihoods, data is 1D or 2D, where axis 0 represents different samples and axis 1 the values. For the Poissonian likelihoods, data is reshaped to align with model_parameters.
model_parameters – Theoretical model values. For the Gaussian likelihoods, this is a 1D array of parameters used as the mean of the multivariate normal. For the Poissonian likelihoods, this is the expected counts (Poisson means).

gaussian(cov: ArrayLike, *, return_log: bool = True) → tuple[tuple[NDArray[float64], ...], NDArray[float64]]#

Evaluates a Gaussian likelihoods for the stored data and parameters.

Parameters:

cov – Covariance matrix. May be a scalar, a 1D vector of diagonal variances, or a full 2D covariance matrix. It will be symmetrized and normalized internally.
return_log – If True, return the log-likelihoods instead of the probability density function.

Returns:

coordinate_grids is a tuple of 1D arrays giving the evaluation coordinates for each dimension.
probabilities is an array with the values of the multivariate Gaussian probability density (or log-density) evaluated on the Cartesian product of those coordinates.

Return type:

A tuple (coordinate_grids, probabilities) where

poissonian(*, return_log: bool = True) → tuple[ndarray, ndarray]#

Evaluates a Poissonian likelihoods for the stored data and parameters.

Parameters:

return_log – If True, return the log-likelihoods instead of the probability mass function.

Returns:

counts is the data reshaped to align with the model parameters.
probabilities is an array of Poisson probabilities (or log-probabilities) computed from counts and model_parameters.

Return type:

A tuple (counts, probabilities) where