bayesreconpy.linear.conditioning#

Conditioning and sampling-based linear reconciliation methods.

bayesreconpy.linear.conditioning.reconc_buis(A, base_forecasts, in_type, distr, num_samples=20000, suppress_warnings=False, seed=None)#

BUIS for Probabilistic Reconciliation of Forecasts via Conditioning

Uses the Bottom-Up Importance Sampling (BUIS) algorithm to draw samples from the reconciled forecast distribution obtained via conditioning.

Parameters:

  • A (ndarray) – Aggregation matrix with shape (n_upper, n_bottom). Each column represents a bottom-level forecast, and each row represents an upper-level forecast. A value of 1 indicates a contribution of a bottom-level forecast to an upper-level forecast.

  • base_forecasts (List[Union[Dict[str, float], ndarray]]) –

    A list containing n_upper + n_bottom elements. The first n_upper elements represent the upper-level base forecasts (in the order of rows in A), and the remaining elements represent the bottom-level base forecasts (in the order of columns in A).

    • If in_type[i] == “samples”, base_forecasts[i] is a NumPy array containing samples from the base forecast distribution.

    • If in_type[i] == “params”, base_forecasts[i] is a dictionary containing parameters for the specified distribution in distr[i]:

      • ’gaussian’: {“mean”: float, “sd”: float}

      • ’poisson’: {“lambda”: float}

      • ’nbinom’: {“size”: float, “prob”: float} or {“size”: float, “mu”: float}

  • in_type (Union[str, List[str]]) –

    Specifies the input type for each base forecast. If a string, the same input type is applied to all forecasts. If a list, in_type[i] specifies the type for the i-th forecast:

    • ’samples’: The forecast is provided as samples.

    • ’params’: The forecast is provided as parameters.

  • distr (Union[str, List[str]]) –

    Specifies the distribution type for each base forecast. If a string, the same distribution is applied to all forecasts. If a list, distr[i] specifies the distribution for the i-th forecast:

    • ’continuous’ or ‘discrete’ if in_type[i] == “samples”.

    • ’gaussian’, ‘poisson’, or ‘nbinom’ if in_type[i] == “params”.

  • num_samples (int) – Number of samples to draw from the reconciled distribution. Ignored if in_type == “samples”, in which case the number of reconciled samples matches the number of samples in the base forecasts. Default is 20,000.

  • suppress_warnings (bool) – Whether to suppress warnings during the importance sampling step. Default is False.

  • seed (Optional[int]) – Random seed for reproducibility. Default is None.

Returns:

A dictionary containing the reconciled forecasts:

  • ’bottom_reconciled_samples’: numpy.ndarray

    A matrix of shape (n_bottom, num_samples) containing the reconciled samples for the bottom-level forecasts.

  • ’upper_reconciled_samples’: numpy.ndarray

    A matrix of shape (n_upper, num_samples) containing the reconciled samples for the upper-level forecasts.

  • ’reconciled_samples’: numpy.ndarray

    A matrix of shape (n, num_samples) containing the reconciled samples for all forecasts.

Return type:

Dict[str, ndarray]

Notes

  • Warnings are triggered during the importance sampling step if:

    • All weights are zero (the corresponding upper forecast is ignored).

    • Effective sample size is less than 200.

    • Effective sample size is less than 1% of the total number of samples.

  • Such warnings indicate potential issues with the base forecasts. Check the inputs in case of warnings.

Examples

Example 1: Gaussian base forecasts (in params form)
>>> A = np.array([
...     [1, 0, 0],
...     [0, 1, 1]
... ])
>>> base_forecasts = [
...     {"mean": 9.0, "sd": 3.0},  # Upper forecast
...     {"mean": 2.0, "sd": 2.0},  # Bottom forecast 1
...     {"mean": 4.0, "sd": 2.0}   # Bottom forecast 2
... ]
>>> result = reconc_buis(A, base_forecasts, in_type="params", distr="gaussian", num_samples=10000, seed=42)
>>> print(result['reconciled_samples'].shape)
(3, 10000)
Example 2: Poisson base forecasts (in params form)
>>> A = np.array([
...     [1, 0, 0],
...     [0, 1, 1]
... ])
>>> base_forecasts = [
...     {"lambda": 9.0},  # Upper forecast
...     {"lambda": 2.0},  # Bottom forecast 1
...     {"lambda": 4.0}   # Bottom forecast 2
... ]
>>> result = reconc_buis(A, base_forecasts, in_type="params", distr="poisson", num_samples=10000, seed=42)
>>> print(result['reconciled_samples'].shape)
(3, 10000)
bayesreconpy.linear.conditioning.reconc_mcmc(A, base_forecasts, distr, num_samples=10000, tuning_int=100, init_scale=1, burn_in=1000, seed=None)#

MCMC for Probabilistic Reconciliation of forecasts via conditioning.

Uses a Markov Chain Monte Carlo (MCMC) algorithm to draw samples from the reconciled forecast distribution, obtained via conditioning. This implementation uses the Metropolis-Hastings algorithm.

Note: This implementation only supports Poisson or Negative Binomial base forecasts. Gaussian distributions are not supported for MCMC.

Parameters:

  • A (ndarray) – Aggregation matrix of shape (n_upper, n_bottom). Each column represents a bottom-level forecast, and each row represents an upper-level forecast. A value of 1 in A[i, j] indicates that bottom-level forecast j contributes to upper-level forecast i.

  • base_forecasts (list) –

    A list containing the parameters of the base forecast distributions. The first n_upper elements correspond to the upper base forecasts (in the order of rows in A), and the remaining n_bottom elements correspond to the bottom base forecasts (in the order of columns in A). Each element is a dictionary:

    • For ‘poisson’: {“lambda”: float}

    • For ‘nbinom’: {“size”: float, “prob”: float} or {“size”: float, “mu”: float}

  • distr (str) –

    Type of predictive distribution. Supported values:

    • ’poisson’

    • ’nbinom’

  • num_samples (int) – Number of samples to draw using MCMC. Default is 10,000.

  • tuning_int (int) – Number of iterations between scale updates of the proposal. Default is 100.

  • init_scale (float) – Initial scale of the proposal distribution. Default is 1.0.

  • burn_in (int) – Number of initial samples to discard. Default is 1,000.

  • seed (Optional[int]) – Random seed for reproducibility. Default is None.

Returns:

A dictionary containing the reconciled forecasts:

  • ’bottom_reconciled_samples’: numpy.ndarray

    A matrix of shape (n_bottom, num_samples) containing reconciled samples for the bottom-level time series.

  • ’upper_reconciled_samples’: numpy.ndarray

    A matrix of shape (n_upper, num_samples) containing reconciled samples for the upper-level time series.

  • ’reconciled_samples’: numpy.ndarray

    A matrix of shape (n, num_samples) containing the reconciled samples for all time series, where n = n_upper + n_bottom.

Return type:

Dict[str, ndarray]

Notes

  • This is a bare-bones implementation of the Metropolis-Hastings algorithm.

  • We recommend using additional tools to assess the convergence of the MCMC chains.

  • The reconc_buis function is generally faster for most hierarchies.

Examples

Example: Simple hierarchy with Poisson base forecasts
>>> import numpy as np
>>> from bayesreconpy.reconc_buis import reconc_buis
>>>
>>> # Create a minimal hierarchy with 1 upper and 2 bottom variables
>>> A = np.array([[1, 1]])  # Aggregation matrix
>>>
>>> # Set the parameters of the Poisson base forecast distributions
>>> lambda_vals = [9, 2, 4]
>>> base_forecasts = [{"lambda": lam} for lam in lambda_vals]
>>>
>>> # Perform MCMC reconciliation
>>> mcmc_result = reconc_mcmc(A, base_forecasts, distr="poisson", num_samples=30000, seed=42)
>>>
>>> # Access reconciled samples
>>> samples_mcmc = mcmc_result['reconciled_samples']
>>>
>>> # Compare reconciled means with those from BUIS
>>> buis_result = reconc_buis(A, base_forecasts, in_type="params", distr="poisson", num_samples=100000, seed=42)
>>> samples_buis = buis_result['reconciled_samples']
>>>
>>> print("MCMC Reconciled Means:", np.mean(samples_mcmc, axis=1))
>>> print("BUIS Reconciled Means:", np.mean(samples_buis, axis=1))

References

  • Corani, G., Azzimonti, D., Rubattu, N. (2024). Probabilistic reconciliation of count time series. International Journal of Forecasting, 40(2), 457-469.

See also

reconc_buis

Faster reconciliation method for most hierarchies.

bayesreconpy.linear.conditioning.reconc_mix_cond(A, fc_bottom, fc_upper, bottom_in_type='pmf', distr=None, num_samples=20000, return_type='pmf', suppress_warnings=False, seed=None)#

Probabilistic forecast reconciliation of mixed hierarchies via conditioning.

Uses importance sampling to draw samples from the reconciled forecast distribution, obtained via conditioning, in the case of a mixed hierarchy.

Parameters:

  • A (ndarray) – Aggregation matrix with shape (n_upper, n_bottom). Each column represents a bottom-level forecast, and each row represents an upper-level forecast. A value of 1 in A[i, j] indicates that bottom-level forecast j contributes to upper-level forecast i.

  • fc_bottom (Union[array, dict]) –

    Base bottom forecasts. The format depends on bottom_in_type:

    • If bottom_in_type == “pmf”: A dictionary where each value is a PMF object (probability mass function).

    • If bottom_in_type == “samples”: A dictionary or array where each value contains samples.

    • If bottom_in_type == “params”: A dictionary where each value contains parameters:

      • ’poisson’: {“lambda”: float}

      • ’nbinom’: {“size”: float, “prob”: float} or {“size”: float, “mu”: float}

  • fc_upper (dict) –

    Base upper forecasts, represented as parameters of a multivariate Gaussian distribution:

    • ’mu’: A vector of length n_upper containing the means.

    • ’Sigma’: A covariance matrix of shape (n_upper, n_upper).

  • bottom_in_type (str) –

    Specifies the type of the base bottom forecasts. Possible values are:

    • ’pmf’: Bottom base forecasts are provided as PMF objects.

    • ’samples’: Bottom base forecasts are provided as samples.

    • ’params’: Bottom base forecasts are provided as estimated parameters.

    Default is ‘pmf’.

  • distr (Optional[str]) –

    Specifies the distribution type for the bottom base forecasts. Only used if bottom_in_type == “params”. Possible values:

    • ’poisson’

    • ’nbinom’

  • num_samples (int) – Number of samples to draw from the reconciled distribution. Ignored if bottom_in_type == “samples”. In this case, the number of reconciled samples equals the number of samples in the base forecasts. Default is 20,000.

  • return_type (str) –

    Specifies the return type of the reconciled distributions. Possible values:

    • ’pmf’: Returns reconciled marginal PMF objects.

    • ’samples’: Returns reconciled multivariate samples.

    • ’all’: Returns both PMF objects and samples.

    Default is ‘pmf’.

  • suppress_warnings (bool) – Whether to suppress warnings about importance sampling weights. Default is False.

  • seed (Optional[int]) – Random seed for reproducibility. Default is None.

Returns:

A dictionary containing the reconciled forecasts:

  • ’bottom_reconciled’: Contains the reconciled forecasts for the bottom-level variables.

    • If return_type == “pmf”: A list of PMF objects.

    • If return_type == “samples”: A matrix of shape (n_bottom, num_samples).

    • If return_type == “all”: Contains both PMF objects and samples.

  • ’upper_reconciled’: Contains the reconciled forecasts for the upper-level variables.

    • If return_type == “pmf”: A list of PMF objects.

    • If return_type == “samples”: A matrix of shape (n_upper, num_samples).

    • If return_type == “all”: Contains both PMF objects and samples.

  • ’ESS’: Effective sample size after importance sampling.

Return type:

Dict[str, Union[Dict[str, Union[List, ndarray]], float]]

Notes

  • A PMF object is a numerical vector where each element corresponds to the probability of integers from 0 to the last value in the support.

  • Warnings are triggered during the importance sampling step if:

    • All weights are zero, causing the upper forecast to be ignored.

    • Effective sample size (ESS) is less than 200.

    • ESS is less than 1% of the total sample size.

Examples

Simple hierarchy with PMF inputs
>>> import numpy as np
>>> from scipy.stats import poisson
>>>
>>> # Simple hierarchy with one upper and two bottom nodes
>>> A = np.array([[1, 1]])  # Aggregation matrix
>>>
>>> # Bottom forecasts as Poisson distributions
>>> lambda_val = 15
>>> n_tot = 60
>>> fc_bottom = {
...     0: np.array([poisson.pmf(x, lambda_val) for x in range(n_tot + 1)]),
...     1: np.array([poisson.pmf(x, lambda_val) for x in range(n_tot + 1)])
... }
>>>
>>> # Upper forecast as Gaussian parameters
>>> fc_upper = {
...     "mu": np.array([40]),  # Mean
...     "Sigma": np.array([[25]])  # Variance (standard deviation squared)
... }
>>>
>>> # Perform reconciliation
>>> result = reconc_mix_cond(A, fc_bottom, fc_upper, bottom_in_type="pmf", return_type="all")
>>>
>>> # Check the results
>>> print(result['bottom_reconciled']['pmf'][0])
>>> print(result['bottom_reconciled']['pmf'][1])
>>> print(result['upper_reconciled']['pmf'][0])

References

  • Zambon, L., Azzimonti, D., Rubattu, N., Corani, G. (2024). Probabilistic reconciliation of mixed-type hierarchical time series. 40th Conference on Uncertainty in Artificial Intelligence.

See also

reconc_td_cond

Probabilistic reconciliation via top-down conditioning.

reconc_buis

Reconciliation using Bottom-Up Importance Sampling (BUIS).

bayesreconpy.linear.conditioning.reconc_td_cond(A, fc_bottom, fc_upper, bottom_in_type='pmf', distr=None, num_samples=20000, return_type='pmf', suppress_warnings=False, seed=None, return_num_samples_ok=False)#

Probabilistic forecast reconciliation of mixed hierarchies via top-down conditioning.

Uses the top-down conditioning algorithm to draw samples from the reconciled forecast distribution. Reconciliation is performed in two steps: 1. Upper base forecasts are reconciled via conditioning using the hierarchical constraints between upper variables. 2. Bottom distributions are updated via a probabilistic top-down procedure.

Parameters:

  • A (ndarray) – Aggregation matrix with shape (n_upper, n_bottom). Each column represents a bottom-level forecast, and each row represents an upper-level forecast. A value of 1 in A[i, j] indicates that bottom-level forecast j contributes to upper-level forecast i.

  • fc_bottom (Union[array, dict]) –

    Base bottom forecasts. The format depends on bottom_in_type:

    • If bottom_in_type == “pmf”: A dictionary where each value is a PMF object (probability mass function).

    • If bottom_in_type == “samples”: A dictionary or array where each value contains samples.

    • If bottom_in_type == “params”: A dictionary where each value contains parameters:

      • ’poisson’: {“lambda”: float}

      • ’nbinom’: {“size”: float, “prob”: float} or {“size”: float, “mu”: float}

  • fc_upper (dict) –

    Base upper forecasts, represented as parameters of a multivariate Gaussian distribution:

    • ’mu’: A vector of length n_upper containing the means.

    • ’Sigma’: A covariance matrix of shape (n_upper, n_upper).

  • bottom_in_type (str) –

    Specifies the type of the base bottom forecasts. Possible values are:

    • ’pmf’: Bottom base forecasts are provided as PMF objects.

    • ’samples’: Bottom base forecasts are provided as samples.

    • ’params’: Bottom base forecasts are provided as estimated parameters.

    Default is ‘pmf’.

  • distr (Optional[str]) –

    Specifies the distribution type for the bottom base forecasts. Only used if bottom_in_type == “params”. Possible values:

    • ’poisson’

    • ’nbinom’

  • num_samples (int) – Number of samples to draw from the reconciled distribution. Ignored if bottom_in_type == “samples”. In this case, the number of reconciled samples equals the number of samples in the base forecasts. Default is 20,000.

  • return_type (str) –

    Specifies the return type of the reconciled distributions. Possible values:

    • ’pmf’: Returns reconciled marginal PMF objects.

    • ’samples’: Returns reconciled multivariate samples.

    • ’all’: Returns both PMF objects and samples.

    Default is ‘pmf’.

  • suppress_warnings (bool) – Whether to suppress warnings about samples lying outside the support of the bottom-up distribution. Default is False.

  • seed (Optional[int]) – Random seed for reproducibility. Default is None.

  • return_num_samples_ok (bool) – Whether to return the number of samples that were valid after reconciliation. Default is False.

Returns:

A dictionary containing the reconciled forecasts:

  • ’bottom_reconciled’: Contains the reconciled forecasts for the bottom-level variables.

    • If return_type == “pmf”: A list of PMF objects.

    • If return_type == “samples”: A matrix of shape (n_bottom, num_samples).

    • If return_type == “all”: Contains both PMF objects and samples.

  • ’upper_reconciled’: Contains the reconciled forecasts for the upper-level variables.

    • If return_type == “pmf”: A list of PMF objects.

    • If return_type == “samples”: A matrix of shape (n_upper, num_samples).

    • If return_type == “all”: Contains both PMF objects and samples.

  • If return_num_samples_ok is True, a tuple is returned.

  • The first element is the above dictionary.

  • The second element is the number of valid samples after reconciliation.

Return type:

Union[Dict[str, Dict[str, Union[List, ndarray]]], Tuple[Dict, int]]

Notes

  • A PMF object is a numerical vector where each element corresponds to the probability of integers from 0 to the last value in the support.

  • Samples lying outside the support of the bottom-up distribution are discarded, and a warning is issued if suppress_warnings is False.

Examples

Simple hierarchy with PMF inputs
>>> import numpy as np
>>> from scipy.stats import poisson
>>>
>>> # Simple hierarchy with one upper and two bottom nodes
>>> A = np.array([[1, 1]])  # Aggregation matrix
>>>
>>> # Bottom forecasts as Poisson distributions
>>> lambda_val = 15
>>> n_tot = 60
>>> fc_bottom = {
...     0: np.array([poisson.pmf(x, lambda_val) for x in range(n_tot + 1)]),
...     1: np.array([poisson.pmf(x, lambda_val) for x in range(n_tot + 1)])
... }
>>>
>>> # Upper forecast as Gaussian parameters
>>> fc_upper = {
...     "mu": np.array([40]),  # Mean
...     "Sigma": np.array([[25]])  # Variance (standard deviation squared)
... }
>>>
>>> # Perform reconciliation
>>> result = reconc_td_cond(A, fc_bottom, fc_upper, bottom_in_type="pmf", return_type="all")
>>>
>>> # Check the results
>>> print(result['bottom_reconciled']['pmf'][0])
>>> print(result['bottom_reconciled']['pmf'][1])
>>> print(result['upper_reconciled']['pmf'][0])

References

  • Zambon, L., Azzimonti, D., Rubattu, N., Corani, G. (2024). Probabilistic reconciliation of mixed-type hierarchical time series. 40th Conference on Uncertainty in Artificial Intelligence.

See also

reconc_mix_cond

Reconciliation for mixed-type hierarchies.

reconc_buis

Reconciliation using Bottom-Up Importance Sampling (BUIS).

bayesreconpy.linear.gaussian#

Gaussian linear reconciliation.

bayesreconpy.linear.gaussian.reconc_gaussian(A, base_forecasts_mu, base_forecasts_Sigma)#

Analytical reconciliation of Gaussian base forecasts. Reconciles forecasts using a closed-form computation for Gaussian base forecasts.

Parameters:

  • A (ndarray) – Aggregation matrix with shape (n_upper, n_bottom). Each column represents a bottom-level forecast, and each row represents an upper-level forecast. A value of 1 in A[i, j] indicates that the bottom-level forecast j contributes to the upper-level forecast i.

  • base_forecasts_mu (List[float]) –

    A 1D list containing the means of the base forecasts. The order is:

    • First, the upper-level means (in the order of rows in A).

    • Then, the bottom-level means (in the order of columns in A).

  • base_forecasts_Sigma (ndarray) – A 2D covariance matrix representing the uncertainties in the base forecasts. The order of rows and columns must match the order of base_forecasts_mu. Shape: (n, n), where n = n_upper + n_bottom.

Returns:

A dictionary containing:

  • ’bottom_reconciled_mean’numpy.ndarray

    A 1D array of the reconciled means for the bottom-level forecasts. Shape: (n_bottom,).

  • ’bottom_reconciled_covariance’numpy.ndarray

    A 2D covariance matrix of the reconciled bottom-level forecasts. Shape: (n_bottom, n_bottom).

Return type:

Dict[str, ndarray]

Notes

  • The function assumes that base forecasts follow a multivariate Gaussian distribution.

  • The covariance matrix base_forecasts_Sigma should be symmetric and positive semi-definite.

  • The order of elements in base_forecasts_mu and rows/columns in base_forecasts_Sigma is critical: first the upper-level forecasts (in the order of rows in A), followed by the bottom-level forecasts (in the order of columns in A).

  • The function returns only the reconciled parameters for the bottom-level forecasts. Reconciled parameters for upper-level forecasts and the entire hierarchy can be derived using the reconciliation matrix A.

Examples

Example 1: Minimal hierarchy with Gaussian base forecasts
>>> A = np.array([
...     [1, 0, 0],
...     [0, 1, 1]
... ])
>>> base_forecasts_mu = [9.0, 2.0, 4.0]
>>> base_forecasts_Sigma = np.diag([9.0, 4.0, 4.0])
>>> result = reconc_gaussian(A, base_forecasts_mu, base_forecasts_Sigma)
>>> print(result['bottom_reconciled_mean'])
[2.5, 4.0]
>>> print(result['bottom_reconciled_covariance'])
[[2.25, 1.5 ],
 [1.5 , 2.0 ]]

References

  • Corani, G., Azzimonti, D., Augusto, J.P.S.C., Zaffalon, M. (2021). Probabilistic Reconciliation of Hierarchical Forecast via Bayes’ Rule. ECML PKDD 2020. Lecture Notes in Computer Science, vol 12459. https://doi.org/10.1007/978-3-030-67664-3_13

  • Zambon, L., Agosto, A., Giudici, P., Corani, G. (2024). Properties of the reconciled distributions for Gaussian and count forecasts. International Journal of Forecasting (in press). https://doi.org/10.1016/j.ijforecast.2023.12.004

bayesreconpy.linear.projection#

Projection-based linear reconciliation methods.

bayesreconpy.linear.projection.estimate_cov_matrix(res, upper=False, n_u=None)#

Estimate the covariance matrix using the Schafer-Strimmer method.

Parameters: - res: array-like, residuals or data for covariance estimation. - upper: bool, whether to use only the upper portion of the data. - n_u: int, number of upper rows/columns to include (required if upper=True).

Returns: - ss: Covariance matrix (shrinkage-based or regular).

bayesreconpy.linear.projection.get_S_from_A(A)#
bayesreconpy.linear.projection.reconc_mint(A, base_forecasts, res=None, W=None, samples=False)#

Reconciles base forecasts using the MinT (Minimum Trace) approach with a shrinkage-based covariance estimator.

The MinT method adjusts base forecasts to ensure coherence with a given aggregation structure, minimizing the total variance of the reconciliation errors using the residual covariance matrix.

Parameters:#

Anp.ndarray

The aggregation constraint matrix (typically of shape [n_agg_levels, n_total_series]), defining how the bottom-level time series aggregate into upper levels.

base_forecastsnp.ndarray

The base forecasts to be reconciled: - of shape [n_total_series, n_time, n_samples] if samples=True, or - of shape [n_total_series, n_time] if samples=False.

resnp.ndarray

Residuals (forecast errors) from a previous model, used to estimate the covariance matrix. Should be of shape [n_total_series, n_time].

Wnp.ndarray, optional

Pre-computed covariance/weight matrix (W_h). If provided, ‘res’ is ignored.

samplesbool, optional (default=False)

Indicates whether the base forecasts include multiple samples. If True, reconciliation is applied sample-by-sample along the third axis.

Returns:#

:

y_tilde_meannp.ndarray

The reconciled forecast means, with the same shape as base_forecasts.

y_tilde_varnp.ndarray or list of np.ndarray

The reconciled forecast variances: - A single matrix of shape [n_total_series, n_total_series] if samples=False. - A list of such matrices, one per sample, if samples=True.

bayesreconpy.linear.projection.reconc_ols(A, base_forecasts, samples=False)#

Reconciles base forecasts using Ordinary Least Squares (OLS) based on a summing matrix.

This function applies forecast reconciliation using the OLS method, ensuring that the reconciled forecasts are coherent with the aggregation structure defined by the matrix A.

Parameters:#

Anp.ndarray

The aggregation constraint matrix (usually of shape [n_agg_levels, n_total_series]). It defines how bottom-level time series are aggregated into upper levels.

base_forecastsnp.ndarray

The base forecasts before reconciliation. Should be: - of shape [n_total_series, n_time, n_samples] if samples=True, or - of shape [n_total_series, n_time] if samples=False.

samplesbool, optional (default=False)

Indicates whether the base forecasts contain multiple samples (i.e., are probabilistic). If True, reconciliation is applied to each sample individually along the third axis.

Returns:#

:

y_tilde_meannp.ndarray

The reconciled forecasts: - of shape [n_total_series, n_time, n_samples] if samples=True, or - of shape [n_total_series, n_time] if samples=False.

bayesreconpy.nonlinear.conditioning#

bayesreconpy.nonlinear.conditioning.reconc_nl_is(f_u, num_samples=10000, n_free=None, seed=None, joint_mean=None, joint_cov=None, eps_cov=1e-09)#

Nonlinear Probabilistic Reconciliation using Importance Sampling (IS).

Uses importance sampling to reconcile forecasts by conditioning on a nonlinear manifold. The method draws samples from the joint distribution of free and constrained variables, computes importance weights based on the constraint function, and resamples to obtain the reconciled distribution.

Parameters:

  • f_u (Callable[[ndarray], ndarray]) –

    The FTC function mapping free variables to constrained variables.

    • Input: array of shape (N, n_free) - N samples of the n_free free variables.

    • Output: array of shape (n_constrained, N) - constrained variables for each sample.

    The function defines the relationship: constrained = f_u(free).

  • num_samples (int, optional) – Number of samples to draw from the reconciled distribution. Default is 10,000.

  • n_free (int) – Number of free (bottom-level) variables. Required parameter. The joint mean and covariance have dimension n_constrained + n_free.

  • seed (int or None, optional) – Random seed for reproducibility. Default is None.

  • joint_mean (numpy.ndarray) –

    Mean vector of the joint distribution over [constrained, free] variables. Shape: (n_constrained + n_free,).

    The first n_constrained elements are the means of the constrained variables, and the remaining n_free elements are the means of the free variables.

  • joint_cov (numpy.ndarray) – Covariance matrix of the joint distribution over [constrained, free] variables. Shape: (n_constrained + n_free, n_constrained + n_free). Must be symmetric and positive semi-definite.

  • eps_cov (float) – Small regularization value added to the diagonal of covariance matrices to ensure numerical stability. Default is 1e-9.

Returns:

A dictionary containing the reconciled forecasts:

  • ’free_reconciled_samples’: numpy.ndarray

    Reconciled samples for the free variables. Shape: (n_free, num_samples).

  • ’constrained_reconciled_samples’: numpy.ndarray

    Reconciled samples for the constrained variables. Shape: (n_constrained, num_samples).

  • ’reconciled_samples’: numpy.ndarray

    Concatenated reconciled samples for all variables [constrained, free]. Shape: (n_constrained + n_free, num_samples).

Return type:

dict

Raises:

  • ValueError – If joint_mean or joint_cov is None. Both parameters are required.

  • TypeError – If joint_mean is None and accessed before validation.

Notes

  • The function computes weights as: w_j ∝ p(U_pred_j, B_j) / p(B_j) where U_pred_j = f_u(B_j) are the predicted constrained variables.

  • Weights are normalized and used to resample particles via multinomial resampling.

  • The joint distribution should reflect the prior beliefs about both free and constrained variables. Ideally, it captures correlations between hierarchical levels.

  • This method is particularly useful for complex nonlinear relationships where the constraint function cannot be inverted analytically.

Examples

Example: Multiple constraints
>>> n_free = 3
>>> n_constrained = 2
>>> n_total = n_free + n_constrained
>>>
>>> joint_mean = np.ones(n_total) * 2.0
>>> joint_cov = np.eye(n_total) * 0.3
>>>
>>> def f_u(B):
...     # B shape: (N, n_free)
...     # Two constraints: sum of first two frees, and third free
...     return np.vstack([
...         np.sum(B[:, :2], axis=1),    # (N,)
...         B[:, 2]                       # (N,)
...     ])  # Output: (n_constrained, N)
>>>
>>> result = reconc_nl_is(f_u, num_samples=5000, n_free=n_free,
...                        joint_mean=joint_mean, joint_cov=joint_cov, seed=42)
>>> print("Reconciled shape:", result["reconciled_samples"].shape)
Reconciled shape: (5, 5000)
bayesreconpy.nonlinear.conditioning.reconc_nl_ukf(free_base_forecasts, in_type, distr, f_u, constrained_base_forecasts, R, num_samples=10000, seed=None)#

Nonlinear Probabilistic Reconciliation using Unscented Kalman Filter (UKF) Conditioning.

Uses the Unscented Kalman Filter to condition the free forecast distributions on the constrained base forecasts via a user-defined nonlinear constraint function. The method performs probabilistic reconciliation by projecting the joint distribution onto the nonlinear manifold defined by f_u(free) = constrained_base_forecasts.

Parameters:

  • free_base_forecasts (List[Dict[str, Union[float, ndarray]]]) –

    A list of n_free dictionaries, one for each free-level forecast.

    • If in_type[i] == “samples”, free_base_forecasts[i] must contain: * ‘samples’: numpy.ndarray of shape (n_samples,) containing samples from the forecast. * ‘residuals’: numpy.ndarray of shape (n_samples,) containing residuals (samples - mean).

    • If in_type[i] == “params”, free_base_forecasts[i] must be a dictionary containing parameters for the specified distribution in distr[i]: * ‘gaussian’: {“mean”: float, “sd”: float}

  • in_type (List[str]) –

    Specifies the input type for each free forecast. Must have length n_free. Each element is one of:

    • ’samples’: The forecast is provided as samples.

    • ’params’: The forecast is provided as parameters.

    Important: All elements must be the same (either all ‘samples’ or all ‘params’). Mixed input types are not supported and will raise NotImplementedError.

  • distr (List[str]) –

    Specifies the distribution type for each free forecast. Must have length n_free.

    • If in_type[i] == “samples”: distr[i] should be descriptive (e.g., ‘gaussian’, ‘continuous’).

    • If in_type[i] == “params”: distr[i] must be one of: * ‘gaussian’: Gaussian distribution

  • f_u (Callable[[ndarray], ndarray]) –

    The FTC function. Maps free variables to constrained variables via the relationship: constrained = f_u(free).

    The function should: - Accept an array of shape (n_free,) representing the free variables. - Return an array of shape (n_constrained,) representing the constrained variables. - Be differentiable (required by the unscented transform).

  • constrained_base_forecasts (ndarray) – Base forecast for the constrained (upper-level) variables. Shape: (n_constrained,). These are the target values that the reconciled free forecasts should aggregate to (according to the constraint function f_u).

  • R (ndarray) – Measurement noise covariance matrix representing uncertainty in the constrained forecasts. Shape: (n_constrained, n_constrained). Must be symmetric and positive semi-definite.

  • num_samples (int) – Number of samples to draw from the reconciled distribution. Default is 10,000.

  • seed (Optional[int]) – Random seed for reproducibility. Default is None.

Returns:

A dictionary containing the reconciled forecasts:

  • ’free_reconciled_samples’: numpy.ndarray

    Reconciled samples for the free (bottom-level) variables. Shape: (n_free, num_samples). Each column is a sample from the reconciled distribution conditioned on the constrained base forecasts.

Return type:

Dict[str, ndarray]

Raises:

  • NotImplementedError – If in_type contains mixed ‘params’ and ‘samples’ entries.

  • ValueError – If required keys are missing from parameter dictionaries or if sample arrays have incompatible shapes.

  • KeyError – If a required parameter key is missing (e.g., ‘sd’ for Gaussian distribution).

Notes

  • The unscented Kalman filter uses sigma points to approximate the mean and covariance of the nonlinearly transformed distributions, avoiding the need for linearization.

  • All free forecasts must use the same input type. For mixed input types, first convert or reconstruct one type before calling this function.

  • The constraint function f_u should be well-defined for the range of the free forecasts.

  • The measurement noise covariance R should reflect the uncertainty in the constrained base forecasts. Smaller values increase the weight of the constrained forecasts.

Examples

ExampleMultiplicative constraint with samples
>>> rng = np.random.default_rng(42)
>>> n_samples = 500
>>> s0 = rng.normal(2.0, 0.2, n_samples)
>>> s1 = rng.normal(3.0, 0.2, n_samples)
>>>
>>> free_base_forecasts = [
...     {"samples": s0, "residuals": s0 - s0.mean()},
...     {"samples": s1, "residuals": s1 - s1.mean()},
... ]
>>> in_type = ["samples"] * 2
>>> distr = ["gaussian"] * 2
>>>
>>> def multiplicative_constraint(free):
...     # free[0] * free[1] should equal constrained value
...     return free[0] * free[1]
>>>
>>> constrained_base_forecasts = np.array([6.0])
>>> R = 0.1 * np.eye(1)
>>>
>>> result = reconc_nl_ukf(free_base_forecasts, in_type, distr,
...                         multiplicative_constraint, constrained_base_forecasts, R,
...                         num_samples=5000, seed=42)
>>> print("Reconciled shape:", result["free_reconciled_samples"].shape)
Reconciled shape: (2, 5000)
bayesreconpy.nonlinear.conditioning.sample_multivariate_gaussian(mean, cov, n_samples, seed=None)#
Return type:

ndarray

bayesreconpy.nonlinear.conditioning.unscented_transform(mu_x, Sigma_x, f, R, alpha=0.001, beta=2, kappa=1)#

bayesreconpy.nonlinear.projection#

bayesreconpy.nonlinear.projection.reconc_nl_ols(samples, f, W=None, n_iter=10, seed=None)#

Nonlinear Probabilistic Reconciliation via Constrained Projection (OLS).

Projects joint forecast samples onto a nonlinear constraint manifold using iterative constrained optimization. The method minimizes the weighted Euclidean distance from original samples to the constraint manifold f(Z) = 0, ensuring all reconciled samples satisfy the nonlinear constraints exactly.

Parameters:

  • samples (ndarray) – Joint forecast samples to be reconciled. Shape: (n_samples, d) where d is the total number of variables (constrained + free variables combined). Each row represents one sample from the joint distribution.

  • f (Callable[[Array], Array]) –

    Constraint function defining the nonlinear reconciliation manifold.

    • Input: JAX array of shape (d,) representing one sample point.

    • Output: JAX array of shape (k,) where k = n_constraints is the number of constraints. The function defines the constraints as f(Z) = 0.

    The constraints express relationships between variables that reconciled samples must satisfy exactly. For example, hierarchical aggregation constraints: top = middle * b2 becomes top - middle * b2 = 0.

  • W (Optional[ndarray]) –

    Weight matrix for the optimization objective. Shape: (d, d). Default is identity matrix np.eye(d).

    Controls the relative importance of different variables when projecting onto the constraint manifold. Can be used to emphasize certain variables or account for different units of measurement.

    The optimization minimizes: ||W @ (Z_proj - Z_orig)||_2^2 subject to f(Z_proj) = 0.

  • n_iter (int) – Number of Newton iterations for solving the constrained optimization problem. Default is 10. Increase for tighter constraint satisfaction or complex manifolds.

  • seed (Optional[int]) – Random seed for reproducibility. Default is None.

Returns:

A dictionary containing the reconciled forecasts and constraint residuals:

  • ’reconciled_samples’: np.ndarray

    Projected samples satisfying the constraint manifold. Shape: (d, n_samples). Each column is a reconciled sample with f(Z_proj) ≈ 0.

  • ’residuals’: np.ndarray

    Constraint residuals f(Z_proj) for each reconciled sample. Shape: (n_constraints, n_samples). Close to zero for successful reconciliation. Larger residuals indicate that the constraint manifold is difficult to satisfy or the optimization did not converge well.

Return type:

Dict[str, ndarray]

Notes

  • Constraint Satisfaction: Unlike sampling-based methods, reconciled samples exactly (or approximately) satisfy the nonlinear constraints f(Z) = 0. This ensures perfect coherency.

  • Optimization Method: Uses constrained Newton iterations (via JNLR library’s make_solver) to project samples onto the manifold. The method solves a constrained least-squares problem at each iteration.

  • Nonlinearity Support: Handles arbitrary smooth nonlinear constraints (multiplicative, exponential, etc.) without linearization assumptions.

  • Weight Matrix: Use W to handle different variable scales or to weight constraints differently. For example, W = np.diag([1.0, 0.1, 0.1, 0.1, 0.1]) prioritizes getting the first variable correct.

  • Convergence: The n_iter parameter controls the Newton iterations. The method should converge for well-posed problems within 10-20 iterations. Check the residuals to assess convergence.

Raises:

  • ValueError – If constraint function input/output shapes are incompatible with samples.

  • TypeError – If samples are not numpy arrays or constraint function is not callable.

Examples

ExampleMultiple constraints with identity weight matrix
>>> rng = np.random.default_rng(42)
>>> n_samples = 2000
>>>
>>> # Create samples
>>> samples = rng.normal(0, 1, (n_samples, 4))
>>>
>>> # Define two constraints
>>> def f(z):
...     return jnp.array([
...         z[0] + z[1] - z[2],         # constraint 1: z0 + z1 = z2
...         z[1] * z[3] - z[0]          # constraint 2: z1 * z3 = z0
...     ])
>>>
>>> result = reconc_nl_ols(samples, f, n_iter=20, seed=42)
>>> Z_proj = result["reconciled_samples"]
>>> residuals = result["residuals"]
>>>
>>> print("Residuals shape:", residuals.shape)
Residuals shape: (2, 2000)
>>> print("Mean abs residuals per constraint:", np.mean(np.abs(residuals), axis=1))
Mean abs residuals per constraint: [1.2e-04 5.8e-05]

bayesreconpy.core#

Core foundational utilities shared by reconciliation methods.

bayesreconpy.utils#

Split utility modules for validation, sampling, constants, and stats.