banquo ====== .. py:module:: banquo .. autoapi-nested-parse:: Banquo. Submodules ---------- .. toctree:: :maxdepth: 1 /autoapi/banquo/banquo/index /autoapi/banquo/kernels/index /autoapi/banquo/numpyro/index Attributes ---------- .. autoapisummary:: banquo.Device banquo.array Exceptions ---------- .. autoapisummary:: banquo.BanquoError banquo.DataMaxExceedsSupportUpperBoundError banquo.DataMinExceedsSupportLowerBoundError banquo.DataRangeExceedsSupportBoundError banquo.MinMaxNormalizationError Classes ------- .. autoapisummary:: banquo.BetaProtocol banquo.Bernstein banquo.NonparanormalBernstein Functions --------- .. autoapisummary:: banquo.add_intercept_column banquo.bernstein_cdf banquo.bernstein_icdf banquo.bernstein_lpdf banquo.bernstein_pdf banquo.chol2inv banquo.device banquo.diag banquo.divide_ns banquo.extract_minmax_parameters banquo.homographic_ns banquo.kahan_sum banquo.logsumexp banquo.minmax_normalization banquo.multi_normal_cholesky_copula_lpdf banquo.multiply_ns banquo.normalize_covariance banquo.shape_handle_wT banquo.shape_handle_wT_posterior banquo.shape_handle_x banquo.std_ns banquo.squared_fractional_graph_laplacian banquo.flat_index banquo.discrete_stochastic_heat_equation_corr banquo.discrete_whittle_matern_fields_corr Package Contents ---------------- .. py:exception:: BanquoError Bases: :py:obj:`Exception` Base class for exceptions in Banquo. .. py:class:: BetaProtocol Bases: :py:obj:`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. :param a: The first shape parameter (alpha) of the Beta distribution. It is an array to allow for vectorized operations over multiple distributions. :type a: array :param b: The second shape parameter (beta) of the Beta distribution. Similar to `a`, it is an array to allow for vectorized operations over multiple distributions. :type b: array .. rubric:: 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. .. rubric:: 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]) .. py:attribute:: a :type: array .. py:attribute:: b :type: array .. py:method:: lpdf(x) Calculate the log probability density function of the beta distribution. .. py:method:: pdf(x) Calculate the probability density function of the beta distribution. .. py:method:: cdf(x) Calculate the cumulative distribution function of the beta distribution. .. py:method:: icdf(x) Calculate the quantile function of the beta distribution. .. py:exception:: DataMaxExceedsSupportUpperBoundError Bases: :py:obj:`DataRangeExceedsSupportBoundError` Exception for data maximum exceeding support's upper bound. .. py:exception:: DataMinExceedsSupportLowerBoundError Bases: :py:obj:`DataRangeExceedsSupportBoundError` Exception for data minimum exceeding support's lower bound. .. py:exception:: DataRangeExceedsSupportBoundError Bases: :py:obj:`MinMaxNormalizationError` Exception for data range exceeding support's boundary. .. py:attribute:: support :type: array .. py:attribute:: data_range :type: array .. py:attribute:: msg_data :type: str .. py:attribute:: msg_support :type: str .. py:data:: Device Type annotation for device objects. .. py:exception:: MinMaxNormalizationError Bases: :py:obj:`BanquoError` Base class for min-max normalization exceptions in Banquo. .. py:function:: add_intercept_column(x, const = 1) Include intercept column to array `x`. The intercept can be any constant number. :param x: Array to add a constant (intercept) column. :type x: array :param const: constant to be included into array `x`, by default 1. :type const: float | int, optional :returns: For an input `x` with dimensions :math:`(n, d)`, it includes a constant column to `x`, resulting in an array with dimensions :math:`(n, d+1)`. :rtype: array .. py:data:: array Type annotation for array objects. For more information, please refer to `array-api `__. .. py:function:: 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, .. math:: 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 :math:`\alpha = j` and :math:`\beta=k-j+1` for :math:`k` basis functions and :math:`j \in \{1, \ldots k\}`. The Bernstein cdf approximation is given by: .. math:: F(x) \approx \mathbf{w}^\mathrm{T} \mathbf{B}(x), with :math:`\mathbf{B}(x) = (B_{1,k}(x), \ldots, B_{k,k}(x))^\mathrm{T}` and weights :math:`\mathbf{w}` are elements in a k-dimensional simplex. :param beta: A class implementing the protocol that defines the interface for a Beta distribution for each basis function, see :class:`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. :type beta: type[BetaProtocol] :param x: 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 :func:`shape_handle_x`. :type x: array :param w: 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 :func:`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. :type w: array :param keepdims: If True, retains reduced axes as dimensions with size 1, by default False. If False, the dimensions are removed. :type keepdims: bool, optional :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. :rtype: array :raises ValueError: If `x` or `w` don't have exactly 5 dimensions. .. py:function:: bernstein_icdf(beta, x, w, keepdims = False) Compute the icdf for a Bernstein-Dirichlet polynomial model. .. py:function:: 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, .. math:: b(x \mid \alpha, \beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \, x^{\alpha-1}(1-x)^{\beta-1}, with :math:`\alpha = j` and :math:`\beta=k-j+1` for :math:`k` basis functions and :math:`j \in \{1, \ldots k\}`. The Bernstein density approximation is given by: .. math:: p(x) \approx \mathbf{w}^\mathrm{T} \mathbf{b}(x), with :math:`\mathbf{b}(x) = (b_{1,k}(x), \ldots, b_{k,k}(x))^\mathrm{T}` and weights :math:`\mathbf{w}` are elements in a k-dimensional simplex. This function returns :math:`\log(\mathbf{w}^\mathrm{T} \mathbf{b}(x))`. :param beta: A class implementing the protocol that defines the interface for a Beta distribution for each basis function, see :class:`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. :type beta: type[BetaProtocol] :param x: 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 :func:`shape_handle_x`. :type x: array :param w: 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 :func:`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. :type w: array :param keepdims: If True, retains reduced axes as dimensions with size 1, by default False. If False, the dimensions are removed. :type keepdims: bool, optional :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. :rtype: array :raises ValueError: If `x` or `w` don't have exactly 5 dimensions. .. rubric:: 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 :func:`BetaProtocol.lpdf` method computes the log-pdf for the inputs in `x`. - The :func:`logsumexp` function is used to aggregate the weighted log-probabilities across the basis functions, ensuring numerical stability. .. py:function:: 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 :func:`bernstein_lpdf`. Considering the beta density function, .. math:: b(x \mid \alpha, \beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \, x^{\alpha-1}(1-x)^{\beta-1}, with :math:`\alpha = j` and :math:`\beta=k-j+1` for :math:`k` basis functions and :math:`j \in \{1, \ldots k\}`. The Bernstein density approximation is given by: .. math:: p(x) \approx \mathbf{w}^\mathrm{T} \mathbf{b}(x), with :math:`\mathbf{b}(x) = (b_{1,k}(x), \ldots, b_{k,k}(x))^\mathrm{T}` and weights :math:`\mathbf{w}` are elements in a k-dimensional simplex. :param beta: A class implementing the protocol that defines the interface for a Beta distribution for each basis function, see :class:`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. :type beta: type[BetaProtocol] :param x: 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 :func:`shape_handle_x`. :type x: array :param w: 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 :func:`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. :type w: array :param keepdims: If True, retains reduced axes as dimensions with size 1, by default False. If False, the dimensions are removed. :type keepdims: bool, optional :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. :rtype: array :raises ValueError: If `x` or `w` don't have exactly 5 dimensions. .. py:function:: chol2inv(spd_chol) Invert a SPD square matrix from its Choleski decomposition. Given a Choleski decomposition :math:`\Sigma` of a matrix :math:`\Sigma`, i.e. :math:`\Sigma = LL^T`, this function returns the inverse :math:`\Sigma^{-1}`. :param spd_chol: Cholesky factor of the correlation/covariance matrix. :type spd_chol: array :returns: Inverse matrix. :rtype: array .. py:function:: 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. :param x: a array object. :type x: array :returns: A device object. :rtype: Device .. rubric:: Notes - `array_api_compat.device `__. - `JIT mechanics `__. .. py:function:: diag(x) Generate a diagonal matrix from array `x`. :param x: One-dimensional array. :type x: array :returns: Diagonal matrix from input array. :rtype: array .. py:function:: 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 :math:`x_1 / x_2`. The function relies on the formula: .. math:: \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)\,. :param x1: Numerator. :type x1: array :param x2: Denominator. :type x2: array :returns: Quotient. :rtype: array .. py:function:: extract_minmax_parameters(x, support = None) Extract the intercept and slop from `x` for the support :math:`[0, 1]`. These parameters can be applied into the linear transformation, given by, .. math:: y = \frac{-a}{b-a} + \frac{1}{b-a} x, to make the data bounded by :math:`[0, 1]`. Where :math:`a` and :math:`b`. are given by: .. math:: a & = \max\{X_{(1)} - \sqrt{S^2/n}, a'\},\\ b & = \min\{X_{(n)} + \sqrt{S^2/n}, b'\}, with :math:`S^2` representing the sample variance, and :math:`X_{(1)}` and :math:`X_{(n)}` denoting the first and last order statistics, respectively. In this formula :math:`x \in [a', b']` :param x: Elements to be transformed. :type x: array :param support: Two-elements array containing the lower and upper bounds for the elements, by default None. If None, `support` is the unbounded interval :math:`(-\infty, \infty)`. :type support: array | None, optional :returns: Two-elements array containing the intercept and slope for a linear transformation. :rtype: array :raises DataRangeExceedsSupportBoundError: If the data range exceeds support's boundary. :raises DataMinExceedsSupportLowerBoundError: If data minimum exceeds support's lower bound. :raises DataMaxExceedsSupportUpperBoundError: If data maximum exceeds support's upper bound. .. note:: See `Unimodal density estimation using Bernstein polynomials `__. .. py:function:: homographic_ns(x) Numerically stable homographic function. Given an array `x`, this function performs a numerically stable calculation of :math:`1/(1+x)`. The function applies :func:`divide_ns` with :math:`x_1 = 1` and :math:`x_2 = 1 + x`. :param x: Elements to extract the homographic function. :type x: array :returns: homographic function, :math:`1/(1+x)`. :rtype: array .. py:function:: 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 :math:`x_1 - x_2`, it is sufficient to use :func:`kahan_sum` with parameters `x1` and `-x2`. :param x1: First term. :type x1: array :param x2: Second term. :type x2: array :returns: Sum of x1 and x2. :rtype: array .. py:function:: logsumexp(x, axis = None, keepdims = False) Compute the log of the sum of exponentials of input elements over a given axis. :param x: Input array to compute the logsumexp over. :type x: array :param axis: Axis or axes along which to perform the reduction, by default None. If None, the operation is performed over all elements. :type axis: int | None, optional :param keepdims: If True, retains reduced axes as dimensions with size 1, by default False. If False, the dimensions are removed. :type keepdims: bool, optional :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. :rtype: array .. rubric:: 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. .. py:function:: minmax_normalization(x, *, support = None, coeffs = None) Transform `x` to the range :math:`[0, 1]`. Linear transform is applied to `x`, given by, .. math:: y = c_1 + c_2 x, See :func:`extract_minmax_parameters` for more information on how :math:`c_1` and :math:`c_2` can be calculated. :param x: Elements to be transformed. :type x: array :param support: Two-elements array containing the lower and upper bounds for the elements, by default None. If None, `support` is the unbounded interval :math:`(-\infty, \infty)`. :type support: array | None, optional :param coeffs: Two-elements array containing the intercept and slope for a linear transformation, by default None. If None, the both parameters will be calculated by :func:`extract_minmax_parameters`. :type coeffs: array | None, optional :returns: Elements in the the range :math:`[0, 1]`. :rtype: array .. py:function:: multi_normal_cholesky_copula_lpdf(marginal, omega_chol) Compute multivariate normal copula lpdf (Cholesky parameterisation). Considering the copula function :math:`C:[0,1]^d\rightarrow [0,1]` and any :math:`(u_1,\dots,u_d)\in[0,1]^d`, such that :math:`u_i = F_i(X_i) = P(X_i \leq x)` are cumulative distribution functions. The multivariate normal copula is given by :math:`C_\Omega(u) = \Phi_\Omega\left(\Phi^{-1}(u_1),\dots, \Phi^{-1}(u_d) \right)`. It is parameterized by the correlation matrix :math:`\Omega = LL^T`, from which :math:`L` is the Cholesky decomposition. Then, the copula density function is given by .. math:: 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 :math:`\log\left(c_\Omega(u)\right)`. :param marginal: Matrix of outcomes from marginal calculations. In this function, :math:`\text{marginal} = \Phi^{-1}(u)`. :type marginal: array :param omega_chol: Cholesky factor of the correlation matrix. :type omega_chol: array :returns: log density of distribution. :rtype: float .. py:function:: 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 :math:`x_1 \times x_2`. The function relies on the formula: .. math:: 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)\,. :param x1: Factor. :type x1: array :param x2: Factor, by default None. If None, the product will be performed through array `x1`. :type x2: array | None, optional :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. :rtype: array | float .. py:function:: normalize_covariance(cov) Normalize a covariance matrix. Assuming a covariance matrix :math:`\Sigma`, the correlation matrix :math:`\Omega` entries are given by: .. math:: \Omega_{ij} = \Sigma_{ij}/\sqrt{\Sigma_{ii} \Sigma_{jj}}\,. :param cov: SPD covariance matrix :type cov: array :returns: SPD correlation matrix. :rtype: array .. py:function:: 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 :func:`bernstein_lpdf`, :func:`bernstein_pdf`, and :func:`bernstein_cdf`, enabling broadcasting. :param w: 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,)`. :type w: array :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)`. :rtype: array :raises ValueError: If `w` has more than 2 dimensions. .. py:function:: 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 :func:`bernstein_lpdf`, :func:`bernstein_pdf`, and :func:`bernstein_cdf`. :param w: 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)`. :type w: array :param chains: Specifies if `w` includes MCMC chain data, by default False. :type chains: bool, optional :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)`. :rtype: array :raises ValueError: If `w` has an incompatible number of dimensions. .. py:function:: 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 :func:`bernstein_lpdf`, :func:`bernstein_pdf`, and :func:`bernstein_cdf`, leveraging broadcasting and vectorized operations. :param x: 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,)`. :type x: array :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)`. :rtype: array :raises ValueError: If `x` has more than 2 dimensions. .. py:function:: 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). :param x: Elements to extract the standard deviation. :type x: array :param axis: 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. :type axis: int | None, optional :param keepdims: If True, retains reduced axes as dimensions with size 1, by default False. If False, the dimensions are removed. :type keepdims: bool, optional :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. :rtype: array | float .. py:function:: 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` :math:`\lambda` with a shifting factor `kappa` and exponent `alpha`, which must both be positive. The fractional operator eigenvalues are given by, .. math:: \widetilde{\lambda} = \left(\kappa^2 \mathbf{I} + \lambda\right)^{\alpha}. :param eigenvalues: Array of eigenvalues of the graph Laplacian matrix. :type eigenvalues: array :param kappa: Shifting factor applied to the eigenvalues, must be positive. :type kappa: float :param alpha: Exponent controlling the fractional power of the transformed Laplacian, must be positive. It has a linear relation with the smoothness. :type alpha: float :returns: Transformed eigenvalues using the squared fractional Laplacian formula. :rtype: array .. rubric:: Notes - `Non-separable Spatio-temporal Graph Kernels via SPDEs `__. .. py:function:: 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. :param node_i: Index of the first node. :type node_i: int :param node_j: Index of the second node. :type node_j: int :param time_i: Index of the time step for the first node. :type time_i: int :param time_j: Index of the time step for the second node. :type time_j: int :param dim_t: Number of time steps for each node. :type dim_t: int :returns: Flattened indices corresponding to the combined node-time positions. :rtype: tuple[int, int] .. py:function:: 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: .. math:: \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, :math:`\mathbf{W}(\xi)` is the derivative, and :math:`\xi = (\mathbf{s}, t) \in \mathcal{D}`. Given the spatial domain :math:`\mathcal{S} \subseteq \mathbb{R}^d` and the time domain :math:`\mathcal{T} \subseteq \mathbb{R}`, the spatiotemporal domain is :math:`\mathcal{D} = \mathcal{S} \times \mathcal{T}`. With :math:`\tau > 0`, the dispersion parameter is represented as :math:`1/\tau`. The thermal diffusivity of the medium in the process :math:`\mathbf{X}(\xi)` is :math:`\gamma > 0`. Here, the graph Laplacian operator :math:`L` is used in place of the continuous one :math:`\Delta`. This correlation function is given by: .. math:: k\left(\lvert t-t'\rvert\right) = \exp\left(-\mathbf{\Gamma} \lvert t-t'\rvert\right). :func:`flat_index` can be used to provide human-friendly access to the matrix elements. :param t: Time indices. :type t: array :param graph_eigenpair: Eigenvalues and eigenvectors of the graph Laplacian. :type graph_eigenpair: tuple[array, array] :param gamma: Medium's (thermal) diffusivity, must be positive. :type gamma: float :param kappa: Shifting factor applied to the spatial eigenvalues of the graph Laplacian, must be positive. :type kappa: float :param alpha: Exponent controlling the fractional power of the transformed Laplacian, must be positive. It has a linear relation with the smoothness. :type alpha: float :returns: Non-separable spatiotemporal correlation matrix. The spatial domain is discretized as a graph. :rtype: array .. rubric:: Notes - `Non-separable Spatio-temporal Graph Kernels via SPDEs `__. .. py:function:: 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 .. math:: \tau \, \left(\kappa^2 - \Delta_{\mathbf{s}}\right)^{\alpha/2} \mathbf{X}(\mathbf{s}) = \mathbf{W}(\mathbf{s}), where :math:`\mathbf{W}(\mathbf{s})` is Gaussian white noise and :math:`\Delta` is the Laplacian operator on a spatial domain :math:`\Omega`. In the discrete setting, the continuous domain :math:`\Omega` is replaced by a tessellated mesh :math:`\Omega_h` with grid spacing :math:`h`, and the continuous Laplacian is approximated by the graph Laplacian :math:`\mathbf{L}`. Assuming the spectral decomposition .. math:: \mathbf{L} = \mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^{\mathrm{T}}, with eigenvalues :math:`\mathbf{\Lambda}` and eigenvectors :math:`\mathbf{Q}`, the covariance matrix of the associated Gaussian Markov random field (GMRF) is given by .. math:: \mathbf{\Sigma} = \mathbf{Q} \tau^{-2}\left(\kappa^2 \mathbf{I} + \mathbf{\Lambda}\right)^{-\alpha} \mathbf{Q}^{\mathrm{T}}, Removing the precision parameter :math:`\tau` via a diagonal scaling produces the correlation matrix .. math:: \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 :math:`\mathbf{S}` are defined by .. math:: \mathbf{S}_{ii} = \sqrt{\sum_{r=1}^d \mathbf{Q}_{ir}^2 \left( \kappa^2 + \lambda_r \right)^{-\alpha}}. :param graph_eigenpair: A tuple (:math:`\mathbf{\Lambda}`, :math:`\mathbf{Q}`) where :math:`\mathbf{\Lambda}` is an array of eigenvalues and :math:`\mathbf{Q}` is the corresponding array of eigenvectors. :type graph_eigenpair: tuple[array, array] :param kappa: A positive shifting factor applied to the eigenvalues. It adjusts the spatial scale of the field. :type kappa: float :param alpha: 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 :math:`d/2`, where :math:`d` is the number of nodes. :type alpha: float :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. :rtype: array .. rubric:: Notes - `Non-separable Spatio-temporal Graph Kernels via SPDEs `__. .. py:class:: Bernstein(weights, *, validate_args = None) Bases: :py:obj:`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). :param weights: 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. :type weights: array :param validate_args: If True, validates the input parameters. By default, None (no validation is applied). :type validate_args: bool | None, optional :raises ValueError: If `weights` is not at least one-dimensional. .. rubric:: Notes - :func:`Bernstein.log_prob` and :func:`Bernstein.cdf` use :func:`banquo.bernstein_lpdf` and :func:`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. .. py:attribute:: arg_constraints .. py:attribute:: reparametrized_params :value: ['weights'] .. py:attribute:: support .. py:attribute:: weights .. py:method:: 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. :param key: The rng_key to be used for generating random numbers. :type key: random.PRNGKey :param sample_shape: The desired shape of independent samples to draw, by default (). If it is set to an empty tuple, it corresponds to a single sample. :type sample_shape: tuple[int, ...], optional :returns: Samples drawn from the distribution with shape `sample_shape + self.batch_shape + self.event_shape`. :rtype: array :raises ValueError: If the provided key is not a valid PRNG key. .. py:method:: log_prob(value) Compute the lpdf of `value` using the Bernstein polynomial model. :param value: 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. :type value: array :returns: The log-probability density function of `value` under the Bernstein model. :rtype: array .. rubric:: Notes The :func:`banquo.bernstein_lpdf` function is used to calculate the lpdf, based on Beta distribution for Bernstein basis functions. .. py:method:: cdf(value) Compute the cdf of `value` using the Bernstein polynomial model. :param value: 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. :type value: array :returns: The cumulative distribution function of `value` under the Bernstein model. :rtype: array .. rubric:: Notes The :func:`banquo.bernstein_cdf` function is used to calculate the cdf, based on Beta distribution for Bernstein basis functions. .. py:method:: icdf(value) Compute the icdf of `value` using the Bernstein polynomial model. :param value: 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. :type value: array :returns: The quantile function of `value` under the Bernstein model. :rtype: array .. rubric:: Notes The :func:`banquo.bernstein_icdf` function is used to calculate the icdf, based on Beta distribution for Bernstein basis functions. .. py:class:: NonparanormalBernstein(weights, correlation_matrix = None, correlation_cholesky = None, *, validate_args = None) Bases: :py:obj:`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: .. math:: 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 :math:`\mathbf{\Psi}(\mathbf{x}) = (\Psi_1(x_1), \ldots, \Psi_d(x_d))^\mathrm{T}`, where :math:`\Psi_i = \Phi^{-1} \circ F_i`. In this case, :math:`F_i` are Bernstein cdf estimators given by: .. math:: F_i(x) \approx \mathbf{w_i}^\mathrm{T} \mathbf{B}(x), and the marginal density estimators are given by: .. math:: p_i(x) \approx \mathbf{w_i}^\mathrm{T} \mathbf{\beta}(x). Both functions are implemented by :class:`Bernstein` model. :param weights: 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. :type weights: array :param correlation_matrix: Correlation matrix of coupling multivariate normal distribution, by default None. :type correlation_matrix: array | None, optional :param correlation_cholesky: Correlation Cholesky factor of coupling multivariate normal distribution, by default None. :type correlation_cholesky: array | None, optional :param validate_args: If True, validates the input parameters. By default, None (no validation is applied). :type validate_args: bool | None, optional :raises ValueError: If `weights` is not at least one-dimensional. .. py:attribute:: arg_constraints .. py:attribute:: support .. py:attribute:: pytree_data_fields :value: 'weights' .. py:attribute:: weights