banquo

Banquo.

Submodules

• banquo.banquo
• banquo.kernels
• banquo.numpyro

Attributes

Device	Type annotation for device objects.
array	Type annotation for array objects.

Exceptions

BanquoError	Base class for exceptions in Banquo.
DataMaxExceedsSupportUpperBoundError	Exception for data maximum exceeding support's upper bound.
DataMinExceedsSupportLowerBoundError	Exception for data minimum exceeding support's lower bound.
DataRangeExceedsSupportBoundError	Exception for data range exceeding support's boundary.
MinMaxNormalizationError	Base class for min-max normalization exceptions in Banquo.

Classes

BetaProtocol	A protocol that defines the interface for a Beta distribution.
Bernstein	Bernstein polynomial-based model.
NonparanormalBernstein	Nonparanormal model with Bernstein polynomial-based marginals.

Functions

add_intercept_column(x[, const])	Include intercept column to array x.
bernstein_cdf(beta, x, w[, keepdims])	Compute the cdf for a Bernstein-Dirichlet polynomial model.
bernstein_icdf(beta, x, w[, keepdims])	Compute the icdf for a Bernstein-Dirichlet polynomial model.
bernstein_lpdf(beta, x, w[, keepdims])	Compute the lpdf for a Bernstein-Dirichlet polynomial model.
bernstein_pdf(beta, x, w[, keepdims])	Compute the pdf for a Bernstein-Dirichlet polynomial model.
chol2inv(spd_chol)	Invert a SPD square matrix from its Choleski decomposition.
device(x)	Wrap function device from array_api_compat.
diag(x)	Generate a diagonal matrix from array x.
divide_ns(x1, x2)	Numerically stable division.
extract_minmax_parameters(x[, support])	Extract the intercept and slop from x for the support \([0, 1]\).
homographic_ns(x)	Numerically stable homographic function.
kahan_sum(x1, x2)	Element-wise Kahan summation algorithm.
logsumexp(x[, axis, keepdims])	Compute the log of the sum of exponentials of input elements over a given axis.
minmax_normalization(x, *[, support, coeffs])	Transform x to the range \([0, 1]\).
multi_normal_cholesky_copula_lpdf(marginal, omega_chol)	Compute multivariate normal copula lpdf (Cholesky parameterisation).
multiply_ns(x1[, x2])	Numerically stable multiplication.
normalize_covariance(cov)	Normalize a covariance matrix.
shape_handle_wT(w)	Reshape and transpose weights w for use with Bernstein functions.
shape_handle_wT_posterior(w[, chains])	Reshape posterior weights w with optional chain handling for MCMC output.
shape_handle_x(x)	Reshape observations x for use with Bernstein functions.
std_ns(x[, axis, keepdims])	Numerically stable calculation of standard deviation.
squared_fractional_graph_laplacian(eigenvalues, kappa, ...)	Compute the squared fractional graph Laplacian.
flat_index(node_i, node_j, time_i, time_j, dim_t)	Calculate flattened indices for combining node and time indices.
discrete_stochastic_heat_equation_corr(t, ...)	Approximate the discrete stochastic heat equation normalized kernel.
discrete_whittle_matern_fields_corr(graph_eigenpair, ...)	Compute the correlation matrix approximating a Whittle–Matérn field.

Package Contents

exception banquo.BanquoError

Bases: Exception

Base class for exceptions in Banquo.

class banquo.BetaProtocol

Bases: Protocol

A protocol that defines the interface for a Beta distribution.

This protocol outlines the required attributes and methods for working with a Beta distribution, including the log probability density function (lpdf), probability density function (pdf),cumulative distribution function (cdf) and inverse cumulative distribution function (icdf) or quantile function.

Parameters:

• a (array) – The first shape parameter (alpha) of the Beta distribution. It is an array to allow for vectorized operations over multiple distributions.

• b (array) – The second shape parameter (beta) of the Beta distribution. Similar to a, it is an array to allow for vectorized operations over multiple distributions.

Notes

• The Beta distribution is a continuous probability distribution defined on the interval [0, 1].

• The a and b parameters define the shape of the distribution. For instance:

  • a = b = 1 gives a uniform distribution.

  • a > b gives a distribution skewed toward 1.

  • a < b gives a distribution skewed toward 0.

• All methods operate on arrays, allowing for efficient vectorized computation of the log-pdf, pdf, cdf, and icdf across multiple Beta distributions and samples.

Example

>>> import numpy as np
>>> from scipy.stats import beta
>>> class ScipyBeta:
>>>     def __init__(self, a, b):
>>>         self.a = a
>>>         self.b = b
>>>
>>>     def lpdf(self, x):
>>>         return beta(self.a, self.b).logpdf(x)
>>>
>>>     def pdf(self, x):
>>>         return beta(self.a, self.b).pdf(x)
>>>
>>>     def cdf(self, x):
>>>         return beta(self.a, self.b).cdf(x)
>>>
>>>     def icdf(self, x):
>>>         return beta.ppf(self.a, self.b)(x)
>>>
>>> beta_dist = ScipyBeta(a=np.array([2]), b=np.array([5]))
>>> x = np.array([0.2, 0.5])
>>> beta_dist.lpdf(x)
array([ 0.89918526, -0.06453852])

a: array

b: array

lpdf(x)

Calculate the log probability density function of the beta distribution.

Parameters:

x (array)

Return type:

array

pdf(x)

Calculate the probability density function of the beta distribution.

Parameters:

x (array)

Return type:

array

cdf(x)

Calculate the cumulative distribution function of the beta distribution.

Parameters:

x (array)

Return type:

array

icdf(x)

Calculate the quantile function of the beta distribution.

Parameters:

x (array)

Return type:

array

exception banquo.DataMaxExceedsSupportUpperBoundError

Bases: DataRangeExceedsSupportBoundError

Exception for data maximum exceeding support’s upper bound.

exception banquo.DataMinExceedsSupportLowerBoundError

Bases: DataRangeExceedsSupportBoundError

Exception for data minimum exceeding support’s lower bound.

exception banquo.DataRangeExceedsSupportBoundError

Bases: MinMaxNormalizationError

Exception for data range exceeding support’s boundary.

support: array

data_range: array

msg_data: str

msg_support: str

banquo.Device

Type annotation for device objects.

exception banquo.MinMaxNormalizationError

Bases: BanquoError

Base class for min-max normalization exceptions in Banquo.

banquo.add_intercept_column(x, const=1)

Include intercept column to array x.

The intercept can be any constant number.

Parameters:

• x (array) – Array to add a constant (intercept) column.

• const (float | int, optional) – constant to be included into array x, by default 1.

Returns:

For an input x with dimensions \((n, d)\), it includes a constant column to x, resulting in an array with dimensions \((n, d+1)\).

Return type:

array

banquo.array

Type annotation for array objects.

For more information, please refer to array-api.

banquo.bernstein_cdf(beta, x, w, keepdims=False)

Compute the cdf for a Bernstein-Dirichlet polynomial model.

This function evaluates the cdf of a weighted sum of Beta distributions, where each Beta distribution forms a basis function in the Bernstein polynomial. The weights (simplex) for each basis function are specified by w, and the inputs to the Beta distributions are given by x.

Considering the beta cumulative distribution function,

\[B(x \mid \alpha, \beta) = \int_0^x \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \, x^{\alpha-1}(1-x)^{\beta-1} \, \mathrm{d}z,\]

with \(\alpha = j\) and \(\beta=k-j+1\) for \(k\) basis functions and \(j \in \{1, \ldots k\}\). The Bernstein cdf approximation is given by:

\[F(x) \approx \mathbf{w}^\mathrm{T} \mathbf{B}(x),\]

with \(\mathbf{B}(x) = (B_{1,k}(x), \ldots, B_{k,k}(x))^\mathrm{T}\) and weights \(\mathbf{w}\) are elements in a k-dimensional simplex.

Parameters:

• beta (type[BetaProtocol]) –

  A class implementing the protocol that defines the interface for a Beta distribution for each basis function, see BetaProtocol. It should accept two arguments a = j and b = k_j, where:

  • j: index of the basis function.

  • k_j: the complement index for the Beta distribution’s second shape parameter.

• x (array) –

  An array of shape (1, 1, n, 1, d), where:

  • n is the number of samples.

  • d is the number of dimensions.

  Each element represents an observation for fitting the Bernstein polynomial model. The other dimensions assigned with 1 are for array broadcasting, see shape_handle_x().

• w (array) –

  An array of shape (c, s, 1, k, d), where:

  • c is the number of MCMC chains.

  • s is the number of MCMC samples per chain.

  • k is the number of basis functions.

  • d is the number of dimensions.

  The elements of w are the weights assigned to each of the k basis functions. The other dimensions assigned with 1 are for array broadcasting, see shape_handle_wT(). In case of the weights are the product of a MCMC sampling algorithm, i.e., it has dimensions (c, s, 1, k, d), it can be used directly into Bernstein functions.

• keepdims (bool, optional) – If True, retains reduced axes as dimensions with size 1, by default False. If False, the dimensions are removed.

Returns:

The cumulative distribution function evaluated at each observation in x, returned as an array of shape (c, s, n, 1, d), where:

• c is the number of MCMC chains.

• s is the number of MCMC samples per chain.

• n is the number of samples.

• d is the number of dimensions.

Each entry corresponds to the cdf of a sample for a specific dimension in the Bernstein polynomial model with parameters w. In case of keepdims is False, the axes of length one will be removed.

Return type:

array

Raises:

ValueError – If x or w don’t have exactly 5 dimensions.

banquo.bernstein_icdf(beta, x, w, keepdims=False)

Compute the icdf for a Bernstein-Dirichlet polynomial model.

Parameters:

• beta (type[BetaProtocol])

• x (array)

• w (array)

• keepdims (bool)

Return type:

array

banquo.bernstein_lpdf(beta, x, w, keepdims=False)

Compute the lpdf for a Bernstein-Dirichlet polynomial model.

This function evaluates the lpdf of a weighted sum of Beta distributions, where each Beta distribution forms a basis function in the Bernstein polynomial. The weights (simplex) for each basis function are specified by w, and the inputs to the Beta distributions are given by x.

Considering the beta density function,

\[b(x \mid \alpha, \beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \, x^{\alpha-1}(1-x)^{\beta-1},\]

with \(\alpha = j\) and \(\beta=k-j+1\) for \(k\) basis functions and \(j \in \{1, \ldots k\}\). The Bernstein density approximation is given by:

\[p(x) \approx \mathbf{w}^\mathrm{T} \mathbf{b}(x),\]

with \(\mathbf{b}(x) = (b_{1,k}(x), \ldots, b_{k,k}(x))^\mathrm{T}\) and weights \(\mathbf{w}\) are elements in a k-dimensional simplex. This function returns \(\log(\mathbf{w}^\mathrm{T} \mathbf{b}(x))\).

Parameters:

• beta (type[BetaProtocol]) –

  A class implementing the protocol that defines the interface for a Beta distribution for each basis function, see BetaProtocol. It should accept two arguments a = j and b = k_j, where:

  • j: index of the basis function.

  • k_j: the complement index for the Beta distribution’s second shape parameter.

• x (array) –

  An array of shape (1, 1, n, 1, d), where:

  • n is the number of samples.

  • d is the number of dimensions.

  Each element represents an observation for fitting the Bernstein polynomial model. The other dimensions assigned with 1 are for array broadcasting, see shape_handle_x().

• w (array) –

  An array of shape (c, s, 1, k, d), where:

  • c is the number of MCMC chains.

  • s is the number of MCMC samples per chain.

  • k is the number of basis functions.

  • d is the number of dimensions.

  The elements of w are the weights assigned to each of the k basis functions. The other dimensions assigned with 1 are for array broadcasting, see shape_handle_wT(). In case of the weights are the product of a MCMC sampling algorithm, i.e., it has dimensions (c, s, 1, k, d), it can be used directly into Bernstein functions.

• keepdims (bool, optional) – If True, retains reduced axes as dimensions with size 1, by default False. If False, the dimensions are removed.

Returns:

The log-probability density function evaluated at each observation in x, returned as an array of shape (c, s, n, 1, d), where:

• c is the number of MCMC chains.

• s is the number of MCMC samples per chain.

• n is the number of samples.

• d is the number of dimensions.

Each entry corresponds to the lpdf of a sample for a specific dimension in the Bernstein polynomial model with parameters w. In case of keepdims is False, the axes of length one will be removed.

Return type:

array

Raises:

ValueError – If x or w don’t have exactly 5 dimensions.

Notes

• This function leverages broadcasting and reshaping to ensure that the log-pdf of the Beta distributions and the weight vectors can be combined element-wise across dimensions and samples.

• The Beta distributions are parameterized by j and k_j, which vary across the number of basis functions k. The BetaProtocol.lpdf() method computes the log-pdf for the inputs in x.

• The logsumexp() function is used to aggregate the weighted log-probabilities across the basis functions, ensuring numerical stability.

banquo.bernstein_pdf(beta, x, w, keepdims=False)

Compute the pdf for a Bernstein-Dirichlet polynomial model.

This function evaluates the pdf of a weighted sum of Beta distributions, where each Beta distribution forms a basis function in the Bernstein polynomial. The weights (simplex) for each basis function are specified by w, and the inputs to the Beta distributions are given by x. This function exponentiate the bernstein_lpdf().

Considering the beta density function,

\[b(x \mid \alpha, \beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \, x^{\alpha-1}(1-x)^{\beta-1},\]

with \(\alpha = j\) and \(\beta=k-j+1\) for \(k\) basis functions and \(j \in \{1, \ldots k\}\). The Bernstein density approximation is given by:

\[p(x) \approx \mathbf{w}^\mathrm{T} \mathbf{b}(x),\]

with \(\mathbf{b}(x) = (b_{1,k}(x), \ldots, b_{k,k}(x))^\mathrm{T}\) and weights \(\mathbf{w}\) are elements in a k-dimensional simplex.

Parameters:

• beta (type[BetaProtocol]) –

  A class implementing the protocol that defines the interface for a Beta distribution for each basis function, see BetaProtocol. It should accept two arguments a = j and b = k_j, where:

  • j: index of the basis function.

  • k_j: the complement index for the Beta distribution’s second shape parameter.

• x (array) –

  An array of shape (1, 1, n, 1, d), where:

  • n is the number of samples.

  • d is the number of dimensions.

  Each element represents an observation for fitting the Bernstein polynomial model. The other dimensions assigned with 1 are for array broadcasting, see shape_handle_x().

• w (array) –

  An array of shape (c, s, 1, k, d), where:

  • c is the number of MCMC chains.

  • s is the number of MCMC samples per chain.

  • k is the number of basis functions.

  • d is the number of dimensions.

  The elements of w are the weights assigned to each of the k basis functions. The other dimensions assigned with 1 are for array broadcasting, see shape_handle_wT(). In case of the weights are the product of a MCMC sampling algorithm, i.e., it has dimensions (c, s, 1, k, d), it can be used directly into Bernstein functions.

• keepdims (bool, optional) – If True, retains reduced axes as dimensions with size 1, by default False. If False, the dimensions are removed.

Returns:

The probability density function evaluated at each observation in x, returned as an array of shape (c, s, n, 1, d), where:

• c is the number of MCMC chains.

• s is the number of MCMC samples per chain.

• n is the number of samples.

• d is the number of dimensions.

Each entry corresponds to the pdf of a sample for a specific dimension in the Bernstein polynomial model with parameters w. In case of keepdims is False, the axes of length one will be removed.

Return type:

array

Raises:

ValueError – If x or w don’t have exactly 5 dimensions.

banquo.chol2inv(spd_chol)

Invert a SPD square matrix from its Choleski decomposition.

Given a Choleski decomposition \(\Sigma\) of a matrix \(\Sigma\), i.e. \(\Sigma = LL^T\), this function returns the inverse \(\Sigma^{-1}\).

Parameters:

spd_chol (array) – Cholesky factor of the correlation/covariance matrix.

Returns:

Inverse matrix.

Return type:

array

banquo.device(x)

Wrap function device from array_api_compat.

JIT-compiled function include JAX transforms, e.g. DynamicJaxprTracer, that has no device attribute. This function fill this gap by considering this case.

Parameters:

x (array) – a array object.

Returns:

A device object.

Return type:

Device

Notes

• array_api_compat.device.

• JIT mechanics.

banquo.diag(x)

Generate a diagonal matrix from array x.

Parameters:

x (array) – One-dimensional array.

Returns:

Diagonal matrix from input array.

Return type:

array

banquo.divide_ns(x1, x2)

Numerically stable division.

Given two arrays of the same shape, x1 and x2, this function performs a numerically stable element-wise division \(x_1 / x_2\). The function relies on the formula:

\[\frac{x_1}{x_2} = \text{sign}(x_1)\text{sign}(x_2) \exp\left(\log(\lvert x_1\rvert) - \log(\lvert x_2\rvert)\right)\,.\]

Parameters:

• x1 (array) – Numerator.

• x2 (array) – Denominator.

Returns:

Quotient.

Return type:

array

banquo.extract_minmax_parameters(x, support=None)

Extract the intercept and slop from x for the support \([0, 1]\).

These parameters can be applied into the linear transformation, given by,

\[y = \frac{-a}{b-a} + \frac{1}{b-a} x,\]

to make the data bounded by \([0, 1]\). Where \(a\) and \(b\). are given by:

\[\begin{split}a & = \max\{X_{(1)} - \sqrt{S^2/n}, a'\},\\ b & = \min\{X_{(n)} + \sqrt{S^2/n}, b'\},\end{split}\]

with \(S^2\) representing the sample variance, and \(X_{(1)}\) and \(X_{(n)}\) denoting the first and last order statistics, respectively. In this formula \(x \in [a', b']\)

Parameters:

• x (array) – Elements to be transformed.

• support (array | None, optional) – Two-elements array containing the lower and upper bounds for the elements, by default None. If None, support is the unbounded interval \((-\infty, \infty)\).

Returns:

Two-elements array containing the intercept and slope for a linear transformation.

Return type:

array

Raises:

• DataRangeExceedsSupportBoundError – If the data range exceeds support’s boundary.

• DataMinExceedsSupportLowerBoundError – If data minimum exceeds support’s lower bound.

• DataMaxExceedsSupportUpperBoundError – If data maximum exceeds support’s upper bound.

Note

See Unimodal density estimation using Bernstein polynomials.

banquo.homographic_ns(x)

Numerically stable homographic function.

Given an array x, this function performs a numerically stable calculation of \(1/(1+x)\). The function applies divide_ns() with \(x_1 = 1\) and \(x_2 = 1 + x\).

Parameters:

x (array) – Elements to extract the homographic function.

Returns:

homographic function, \(1/(1+x)\).

Return type:

array

banquo.kahan_sum(x1, x2)

Element-wise Kahan summation algorithm.

Given two arrays of the same shape, x1 and x2, this function performs a numerically stable element-wise summation x1 + x2. For the subtraction \(x_1 - x_2\), it is sufficient to use kahan_sum() with parameters x1 and -x2.

Parameters:

• x1 (array) – First term.

• x2 (array) – Second term.

Returns:

Sum of x1 and x2.

Return type:

array

banquo.logsumexp(x, axis=None, keepdims=False)

Compute the log of the sum of exponentials of input elements over a given axis.

Parameters:

• x (array) – Input array to compute the logsumexp over.

• axis (int | None, optional) – Axis or axes along which to perform the reduction, by default None. If None, the operation is performed over all elements.

• keepdims (bool, optional) – If True, retains reduced axes as dimensions with size 1, by default False. If False, the dimensions are removed.

Returns:

An array with the same type as x containing the log of the sum of exponentials. If keepdims is True, the result will have the same number of dimensions as the input. Otherwise, the reduced dimensions are removed.

Return type:

array

Notes

This implementation follows a numerically stable approach to compute the log of the sum of exponentials by:

• Shifting input values by the maximum along the specified axis.

• Computing the exponentials of the shifted values to avoid overflow or underflow.

banquo.minmax_normalization(x, *, support=None, coeffs=None)

Transform x to the range \([0, 1]\).

Linear transform is applied to x, given by,

\[y = c_1 + c_2 x,\]

See extract_minmax_parameters() for more information on how \(c_1\) and \(c_2\) can be calculated.

Parameters:

• x (array) – Elements to be transformed.

• support (array | None, optional) – Two-elements array containing the lower and upper bounds for the elements, by default None. If None, support is the unbounded interval \((-\infty, \infty)\).

• coeffs (array | None, optional) – Two-elements array containing the intercept and slope for a linear transformation, by default None. If None, the both parameters will be calculated by extract_minmax_parameters().

Returns:

Elements in the the range \([0, 1]\).

Return type:

array

banquo.multi_normal_cholesky_copula_lpdf(marginal, omega_chol)

Compute multivariate normal copula lpdf (Cholesky parameterisation).

Considering the copula function \(C:[0,1]^d\rightarrow [0,1]\) and any \((u_1,\dots,u_d)\in[0,1]^d\), such that \(u_i = F_i(X_i) = P(X_i \leq x)\) are cumulative distribution functions. The multivariate normal copula is given by \(C_\Omega(u) = \Phi_\Omega\left(\Phi^{-1}(u_1),\dots, \Phi^{-1}(u_d) \right)\). It is parameterized by the correlation matrix \(\Omega = LL^T\), from which \(L\) is the Cholesky decomposition. Then, the copula density function is given by

\[c_\Omega(u) = \frac{\partial^d C_\Omega(u)}{\partial \Phi(u_1)\cdots \partial \Phi(u_d)} \,,\]

and this function computes its log density \(\log\left(c_\Omega(u)\right)\).

Parameters:

• marginal (array) – Matrix of outcomes from marginal calculations. In this function, \(\text{marginal} = \Phi^{-1}(u)\).

• omega_chol (array) – Cholesky factor of the correlation matrix.

Returns:

log density of distribution.

Return type:

float

banquo.multiply_ns(x1, x2=None)

Numerically stable multiplication.

Given two arrays of the same shape, x1 and x2, this function performs a numerically stable element-wise division \(x_1 \times x_2\). The function relies on the formula:

\[x_1 \times x_2 = \text{sign}(x_1)\text{sign}(x_2) \exp\left(\log(\lvert x_1\rvert) + \log(\lvert x_2\rvert)\right)\,.\]

Parameters:

• x1 (array) – Factor.

• x2 (array | None, optional) – Factor, by default None. If None, the product will be performed through array x1.

Returns:

Product. If x2 is None, the product will be performed in array x1, resulting in a float. If otherwise, the function returns the element-wise multiplication between x1 and x2 resulting in an array of the same shape.

Return type:

array | float

banquo.normalize_covariance(cov)

Normalize a covariance matrix.

Assuming a covariance matrix \(\Sigma\), the correlation matrix \(\Omega\) entries are given by:

\[\Omega_{ij} = \Sigma_{ij}/\sqrt{\Sigma_{ii} \Sigma_{jj}}\,.\]

Parameters:

cov (array) – SPD covariance matrix

Returns:

SPD correlation matrix.

Return type:

array

banquo.shape_handle_wT(w)

Reshape and transpose weights w for use with Bernstein functions.

This function reshapes the weights array w and transposes it for compatibility with functions like bernstein_lpdf(), bernstein_pdf(), and bernstein_cdf(), enabling broadcasting.

Parameters:

w (array) – An array of weights with shape (d, k), where d is the number of dimensions, and k is the number of basis functions. If w is one-dimensional, it represents k basis functions with shape (k,).

Returns:

Reshaped and transposed array of w for Bernstein functions.

• If w has shape (k,), it will be reshaped to (1, 1, 1, k, 1).

• If w has shape (d, k), it will be reshaped to (1, 1, 1, k, d).

Return type:

array

Raises:

ValueError – If w has more than 2 dimensions.

banquo.shape_handle_wT_posterior(w, chains=False)

Reshape posterior weights w with optional chain handling for MCMC output.

This function reshapes and transposes the input weights array w based on whether MCMC chains are included. It ensures compatibility with functions like bernstein_lpdf(), bernstein_pdf(), and bernstein_cdf().

Parameters:

• w (array) – An array of weights with shape (s, k) or (s, d, k) when chains=False. If chains=True, w should have shape (c, s, k) or (c, s, d, k).

• chains (bool, optional) – Specifies if w includes MCMC chain data, by default False.

Returns:

Reshaped array w suitable for Bernstein functions.

• If chains=True, reshapes w to (c, s, 1, k, 1) or (c, s, 1, k, d).

• If chains=False, reshapes w to (1, s, 1, k, 1) or (1, s, 1, k, d).

Return type:

array

Raises:

ValueError – If w has an incompatible number of dimensions.

banquo.shape_handle_x(x)

Reshape observations x for use with Bernstein functions.

This function reshapes the input array x to ensure compatibility with functions such as bernstein_lpdf(), bernstein_pdf(), and bernstein_cdf(), leveraging broadcasting and vectorized operations.

Parameters:

x (array) – An array of observations, with shape (n, d) where n is the number of samples, and d is the number of dimensions. If x is one-dimensional, it represents n samples with shape (n,).

Returns:

Reshaped array of x suitable for Bernstein functions.

• If x has shape (n,), it will be reshaped to (1, 1, n, 1, 1).

• If x has shape (n, d), it will be reshaped to (1, 1, n, 1, d).

Return type:

array

Raises:

ValueError – If x has more than 2 dimensions.

banquo.std_ns(x, axis=None, keepdims=False)

Numerically stable calculation of standard deviation.

If the standard deviation tends to infinity, it is substituted by the interquartile range (IQR).

Parameters:

• x (array) – Elements to extract the standard deviation.

• axis (int | None, optional) – Axis or axes along which the standard deviation calculation is performed. If axis is negative it counts from the last to the first axis, by default None will sum all of the elements of the input array.

• keepdims (bool, optional) – If True, retains reduced axes as dimensions with size 1, by default False. If False, the dimensions are removed.

Returns:

An array with the same shape as a, with the specified axis removed. If a is a 0-d array, or if axis is None, a scalar is returned. If an output array is specified, a reference to out is returned. If keepdims is True, retains reduced axes as dimensions with size 1.

Return type:

array | float

banquo.squared_fractional_graph_laplacian(eigenvalues, kappa, alpha)

Compute the squared fractional graph Laplacian.

This function applies the squared fractional transformation to the graph Laplacian eigenvalues \(\lambda\) with a shifting factor kappa and exponent alpha, which must both be positive. The fractional operator eigenvalues are given by,

\[\widetilde{\lambda} = \left(\kappa^2 \mathbf{I} + \lambda\right)^{\alpha}.\]

Parameters:

• eigenvalues (array) – Array of eigenvalues of the graph Laplacian matrix.

• kappa (float) – Shifting factor applied to the eigenvalues, must be positive.

• alpha (float) – Exponent controlling the fractional power of the transformed Laplacian, must be positive. It has a linear relation with the smoothness.

Returns:

Transformed eigenvalues using the squared fractional Laplacian formula.

Return type:

array

Notes

• Non-separable Spatio-temporal Graph Kernels via SPDEs.

banquo.flat_index(node_i, node_j, time_i, time_j, dim_t)

Calculate flattened indices for combining node and time indices.

Given node indices and time indices, this function returns the equivalent flat indices for a 2D matrix where each node spans dim_t time steps.

Parameters:

• node_i (int) – Index of the first node.

• node_j (int) – Index of the second node.

• time_i (int) – Index of the time step for the first node.

• time_j (int) – Index of the time step for the second node.

• dim_t (int) – Number of time steps for each node.

Returns:

Flattened indices corresponding to the combined node-time positions.

Return type:

tuple[int, int]

banquo.discrete_stochastic_heat_equation_corr(t, graph_eigenpair, gamma, kappa, alpha)

Approximate the discrete stochastic heat equation normalized kernel.

This function builds an approximation for the Gram matrix of the normalized kernel (correlation) resulted from the discrete stochastic heat equation on a graph. The following is the expression for the stochastic heat equation:

\[\left[\frac{\partial}{\partial t} + \gamma \left(\kappa^2 + \Delta\right)^{\alpha/2}\right] \tau \mathbf{X}(\xi) = \mathbf{W}(\xi).\]

For the spatiotemporal Brownian motion, \(\mathbf{W}(\xi)\) is the derivative, and \(\xi = (\mathbf{s}, t) \in \mathcal{D}\). Given the spatial domain \(\mathcal{S} \subseteq \mathbb{R}^d\) and the time domain \(\mathcal{T} \subseteq \mathbb{R}\), the spatiotemporal domain is \(\mathcal{D} = \mathcal{S} \times \mathcal{T}\). With \(\tau > 0\), the dispersion parameter is represented as \(1/\tau\). The thermal diffusivity of the medium in the process \(\mathbf{X}(\xi)\) is \(\gamma > 0\). Here, the graph Laplacian operator \(L\) is used in place of the continuous one \(\Delta\).

This correlation function is given by:

\[k\left(\lvert t-t'\rvert\right) = \exp\left(-\mathbf{\Gamma} \lvert t-t'\rvert\right).\]

flat_index() can be used to provide human-friendly access to the matrix elements.

Parameters:

• t (array) – Time indices.

• graph_eigenpair (tuple[array, array]) – Eigenvalues and eigenvectors of the graph Laplacian.

• gamma (float) – Medium’s (thermal) diffusivity, must be positive.

• kappa (float) – Shifting factor applied to the spatial eigenvalues of the graph Laplacian, must be positive.

• alpha (float) – Exponent controlling the fractional power of the transformed Laplacian, must be positive. It has a linear relation with the smoothness.

Returns:

Non-separable spatiotemporal correlation matrix. The spatial domain is discretized as a graph.

Return type:

array

Notes

• Non-separable Spatio-temporal Graph Kernels via SPDEs.

banquo.discrete_whittle_matern_fields_corr(graph_eigenpair, kappa, alpha)

Compute the correlation matrix approximating a Whittle–Matérn field.

This function constructs an approximation to the kernel corresponding to a Whittle–Matérn field defined as the solution to the SPDE

\[\tau \, \left(\kappa^2 - \Delta_{\mathbf{s}}\right)^{\alpha/2} \mathbf{X}(\mathbf{s}) = \mathbf{W}(\mathbf{s}),\]

where \(\mathbf{W}(\mathbf{s})\) is Gaussian white noise and \(\Delta\) is the Laplacian operator on a spatial domain \(\Omega\). In the discrete setting, the continuous domain \(\Omega\) is replaced by a tessellated mesh \(\Omega_h\) with grid spacing \(h\), and the continuous Laplacian is approximated by the graph Laplacian \(\mathbf{L}\). Assuming the spectral decomposition

\[\mathbf{L} = \mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^{\mathrm{T}},\]

with eigenvalues \(\mathbf{\Lambda}\) and eigenvectors \(\mathbf{Q}\), the covariance matrix of the associated Gaussian Markov random field (GMRF) is given by

\[\mathbf{\Sigma} = \mathbf{Q} \tau^{-2}\left(\kappa^2 \mathbf{I} + \mathbf{\Lambda}\right)^{-\alpha} \mathbf{Q}^{\mathrm{T}},\]

Removing the precision parameter \(\tau\) via a diagonal scaling produces the correlation matrix

\[\mathbf{\Omega} = \mathbf{S}^{-1} \mathbf{Q} \left(\kappa^2 \mathbf{I} + \mathbf{\Lambda}\right)^{-\alpha} \mathbf{Q}^{\mathrm{T}} \mathbf{S}^{-1} ,\]

where the diagonal entries of \(\mathbf{S}\) are defined by

\[\mathbf{S}_{ii} = \sqrt{\sum_{r=1}^d \mathbf{Q}_{ir}^2 \left( \kappa^2 + \lambda_r \right)^{-\alpha}}.\]

Parameters:

• graph_eigenpair (tuple[array, array]) – A tuple (\(\mathbf{\Lambda}\), \(\mathbf{Q}\)) where \(\mathbf{\Lambda}\) is an array of eigenvalues and \(\mathbf{Q}\) is the corresponding array of eigenvectors.

• kappa (float) – A positive shifting factor applied to the eigenvalues. It adjusts the spatial scale of the field.

• alpha (float) – A positive exponent controlling the fractional power of the Laplacian transformation. This parameter is directly related to the smoothness of the field. It must be greater than \(d/2\), where \(d\) is the number of nodes.

Returns:

The normalized correlation matrix approximating the Whittle–Matérn field on the discretized domain. This matrix converges to the true continuous correlation structure as the discretization is refined.

Return type:

array

Notes

• Non-separable Spatio-temporal Graph Kernels via SPDEs.

class banquo.Bernstein(weights, *, validate_args=None)

Bases: numpyro.distributions.Distribution

Bernstein polynomial-based model.

The Bernstein class implements a probability distribution for nonparametric modeling densities with Bernstein polynomials. It uses a set of weights as coefficients for the basis functions and supports computing the log-probability density and cumulative distribution function (CDF).

Parameters:

• weights (array) – Array of weights (simplex) with shape (d, k), where d is the number of dimensions, and k is the number of basis functions (Bernstein polynomial order). If the shape is k, the system will be considered as a one-dimensional array. The weights are, for each dimension d, a k-dimensional unit simplex. The weights are applied as coefficients for the Bernstein polynomial basis in each dimension.

• validate_args (bool | None, optional) – If True, validates the input parameters. By default, None (no validation is applied).

Raises:

ValueError – If weights is not at least one-dimensional.

Notes

• Bernstein.log_prob() and Bernstein.cdf() use banquo.bernstein_lpdf() and banquo.bernstein_cdf() respectively, with Beta distribution based Bernstein basis functions.

• The weights must sum to 1 across each dimension, fulfilling the Dirichlet distribution constraint.

arg_constraints

reparametrized_params = ['weights']

support

weights

sample(key, sample_shape=())

Sample from Bernstein distribution.

Returns a sample from the distribution with the shape specified by sample_shape + batch_shape + event_shape. Exchangeable draws from the distribution instance will fill the leading dimensions (of size sample_shape) of the returned sample when sample_shape is not empty.

Parameters:

• key (random.PRNGKey) – The rng_key to be used for generating random numbers.

• sample_shape (tuple[int, ...], optional) – The desired shape of independent samples to draw, by default (). If it is set to an empty tuple, it corresponds to a single sample.

Returns:

Samples drawn from the distribution with shape sample_shape + self.batch_shape + self.event_shape.

Return type:

array

Raises:

ValueError – If the provided key is not a valid PRNG key.

log_prob(value)

Compute the lpdf of value using the Bernstein polynomial model.

Parameters:

value (array) – The observed data or values to evaluate the lpdf. The array should have shape (n, d), where n is the number of samples, and d is the number of dimensions. If value is one-dimensional with shape (n,), it will be reshaped to (n, 1). Each element represents a sample to be evaluated under the Bernstein polynomial model.

Returns:

The log-probability density function of value under the Bernstein model.

Return type:

array

Notes

The banquo.bernstein_lpdf() function is used to calculate the lpdf, based on Beta distribution for Bernstein basis functions.

cdf(value)

Compute the cdf of value using the Bernstein polynomial model.

Parameters:

value (array) – The observed data or values to evaluate the cdf. The array should have shape (n, d), where n is the number of samples, and d is the number of dimensions. If value is one-dimensional with shape (n,), it will be reshaped to (n, 1). Each element represents a sample to be evaluated under the Bernstein polynomial model.

Returns:

The cumulative distribution function of value under the Bernstein model.

Return type:

array

Notes

The banquo.bernstein_cdf() function is used to calculate the cdf, based on Beta distribution for Bernstein basis functions.

icdf(value)

Compute the icdf of value using the Bernstein polynomial model.

Parameters:

value (array) – The observed data or values to evaluate the icdf. The array should have shape (n, d), where n is the number of samples, and d is the number of dimensions. If value is one-dimensional with shape (n,), it will be reshaped to (n, 1). Each element represents a value in the interval [0, 1] to return quantile under the Bernstein polynomial model.

Returns:

The quantile function of value under the Bernstein model.

Return type:

array

Notes

The banquo.bernstein_icdf() function is used to calculate the icdf, based on Beta distribution for Bernstein basis functions.

class banquo.NonparanormalBernstein(weights, correlation_matrix=None, correlation_cholesky=None, *, validate_args=None)

Bases: numpyro.distributions.GaussianCopula

Nonparanormal model with Bernstein polynomial-based marginals.

The NonparanormalBernstein class implements a probability distribution for nonparametric marginal densities modeling with Bernstein polynomials. It uses a set of weights as coefficients for the basis functions and supports computing the log-probability density and cumulative distribution function (CDF). A Gaussian copula is used to create the joint distribution. The nonparanormal model is formed up of both the nonparametric marginals and the Gaussian copula, and its probability density function is as follows:

\[p(\mathbf{x}) = \lvert\mathbf{\Sigma}\rvert^{-1/2} \, \exp \left( -\frac{1}{2} \mathbf{\Psi}(\mathbf{x})^\mathrm{T} \, \left(\mathbf{\Sigma}^{-1} - \mathbf{I}\right) \, \mathbf{\Psi}(\mathbf{x})\right) \prod_{i=1}^d \lvert p(x_i)\rvert,\]

with \(\mathbf{\Psi}(\mathbf{x}) = (\Psi_1(x_1), \ldots, \Psi_d(x_d))^\mathrm{T}\), where \(\Psi_i = \Phi^{-1} \circ F_i\). In this case, \(F_i\) are Bernstein cdf estimators given by:

\[F_i(x) \approx \mathbf{w_i}^\mathrm{T} \mathbf{B}(x),\]

and the marginal density estimators are given by:

\[p_i(x) \approx \mathbf{w_i}^\mathrm{T} \mathbf{\beta}(x).\]

Both functions are implemented by Bernstein model.

Parameters:

• weights (array) – Array of weights (simplex) with shape (d, k), where d is the number of dimensions, and k is the number of basis functions (Bernstein polynomial order). If the shape is k, the system will be considered as a one-dimensional array. The weights are, for each dimension d, a k-dimensional unit simplex. The weights are applied as coefficients for the Bernstein polynomial basis in each dimension.

• correlation_matrix (array | None, optional) – Correlation matrix of coupling multivariate normal distribution, by default None.

• correlation_cholesky (array | None, optional) – Correlation Cholesky factor of coupling multivariate normal distribution, by default None.

• validate_args (bool | None, optional) – If True, validates the input parameters. By default, None (no validation is applied).

Raises:

ValueError – If weights is not at least one-dimensional.

arg_constraints

support

pytree_data_fields = 'weights'

weights