Dev bugs v1

DESCRIPTION

@@ -57,6 +57,7 @@ Suggests:
     forecast,
     glarma,
     scoringRules,
+    ggplot2,
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
 VignetteBuilder: knitr

R/PMF.R

@@ -141,7 +141,7 @@ PMF_from_params <- function(params, distr, Rtoll = .RTOLL) {
 #'
 #' @param pmf A PMF object (numeric vector where element j+1 is the probability of value j).
 #' @param N_samples Number of samples to draw from the PMF.
-#' @param p Probability level for quantile computation (between 0 and 1).
+#' @param p Probability level for quantile computation. Either a number or a vector with all elements between 0 and 1.
 #' @param Ltoll Lower tolerance for computing the minimum of the PMF (default: 1e-15).
 #' @param Rtoll Upper tolerance for computing the maximum of the PMF (default: 1e-9).
 #'
@@ -189,6 +189,7 @@ PMF_from_params <- function(params, distr, Rtoll = .RTOLL) {
 #'     Computes the quantile of the distribution at probability level \code{p}.
 #'     Returns the smallest value \code{x} such that the cumulative probability up to \code{x}
 #'     is greater than or equal to \code{p}. For example, \code{p=0.5} gives the median.
+#'     The parameter \code{p} can be a single number or a vector of probabilities.
 #'   }
 #'
 #'   \item{\code{PMF_summary(pmf, Ltoll, Rtoll)}}{
@@ -225,11 +226,11 @@ PMF_get_var <- function(pmf) {
 #' @rdname PMF
 #' @export
 PMF_get_quantile <- function(pmf, p) {
-  if (p <= 0 | p >= 1) {
+  if (min(p) <= 0 | max(p) >= 1) {
     stop("Input error: probability p must be between 0 and 1")
   }
   cdf <- cumsum(pmf)
-  x <- min(which(cdf >= p))
+  x <- pmin(findInterval(p, cdf) + 1, length(cdf))
   return(x - 1)
 }
 

R/bayesRecon-package.R

@@ -5,20 +5,20 @@
 #' @section Reconciliation functions:
 #'
 #' The package implements reconciliation via conditioning for probabilistic forecasts
-#' of hierarchical time series. The main functions are:
+#' of hierarchical time series. The reconciliation functions are:
 #'
-#' * [reconc_gaussian()]: reconciliation via conditioning, assumes multivariate Gaussian
+#' * [reconc_gaussian()]: reconciliation via conditioning assuming multivariate Gaussian
 #'    base forecasts; this is done analytically;
+#' * [reconc_t()]: reconciliation via conditioning assuming multivariate Gaussian
+#'    base forecasts with uncertain covariance matrix; 
+#'    the reconciled forecasts are multivariate Student-t; this is done analytically;
 #' * [reconc_BUIS()]: reconciliation via conditioning of any probabilistic forecast
-#'    via importance sampling; this is the recommended option for non-Gaussian base forecasts;
-#' * [reconc_MCMC()]: reconciliation via conditioning of discrete probabilistic
-#'    forecasts via Markov Chain Monte Carlo;
-#' * [reconc_MixCond()]: reconciliation via conditioning of mixed hierarchies, where
+#'    via bottom-up importance sampling; an alternative method for discrete forecasts 
+#'    is implemented in [reconc_MCMC()], but we recommend using `reconc_BUIS`;
+#' * [reconc_MixCond()] and [reconc_TDcond()]: reconciliation of mixed hierarchies, where
 #'    the upper forecasts are multivariate Gaussian and the bottom forecasts are discrete distributions;
-#' * [reconc_TDcond()]: reconciliation via top-down conditioning of mixed hierarchies, where
-#'    the upper forecasts are multivariate Gaussian and the bottom forecasts are discrete distributions;
-#' * [reconc_t()]: reconciliation via conditioning with uncertain covariance matrix; 
-#'    the reconciled forecasts are multivariate Student-t; this is done analytically.
+#'    `reconc_MixCond` implements conditioning via importance sampling, 
+#'    while `reconc_TDcond` implements top-down conditioning.
 #'
 #' @section Utility functions:
 #'

R/reconc_TDcond.R

@@ -287,7 +287,7 @@ reconc_TDcond <- function(A, base_fc_bottom, base_fc_upper,
     U <- .MVN_sample(
       n_samples = num_samples,
       mu = rec_gauss_u$bottom_rec_mean,
-      Sigma = rec_gauss_u$bottom_rec_covariance
+      Sigma = rec_gauss_u$bottom_rec_cov
     )
     U <- round(U) # round
     mode(U) <- "integer" # convert to integer

R/reconc_gaussian.R

@@ -23,9 +23,9 @@
 #' @return A list containing the reconciled forecasts. The list has the following named elements:
 #'
 #' * `bottom_rec_mean`: reconciled mean for the bottom forecasts;
-#' * `bottom_rec_covariance`: reconciled covariance for the bottom forecasts;
+#' * `bottom_rec_cov`: reconciled covariance for the bottom forecasts;
 #' * `upper_rec_mean`: (only if `return_upper = TRUE`) reconciled mean for the upper forecasts;
-#' * `upper_rec_covariance`: (only if `return_upper = TRUE`) reconciled covariance for the upper forecasts.
+#' * `upper_rec_cov`: (only if `return_upper = TRUE`) reconciled covariance for the upper forecasts.
 #'
 #'
 #' @examples
@@ -55,7 +55,7 @@
 #' )
 #'
 #' bottom_mean_rec <- analytic_rec$bottom_rec_mean
-#' bottom_cov_rec <- analytic_rec$bottom_rec_covariance
+#' bottom_cov_rec <- analytic_rec$bottom_rec_cov
 #'
 #' # To obtain reconciled samples for the entire hierarchy, sample from the reconciled 
 #' # bottom distribution and then aggregate using A. 
@@ -110,9 +110,9 @@
 #'   result <- reconc_gaussian(A, base_fc_mean = base_fc_mean, residuals = res, return_upper = TRUE)
 #'
 #'   bottom_mean <- result$bottom_rec_mean
-#'   bottom_cov <- result$bottom_rec_covariance
+#'   bottom_cov <- result$bottom_rec_cov
 #'   upper_mean <- result$upper_rec_mean
-#'   upper_cov <- result$upper_rec_covariance
+#'   upper_cov <- result$upper_rec_cov
 #'
 #'   # Print reconciled means
 #'   cat("Reconciled bottom means:", round(bottom_mean, 3), "\n")
@@ -203,12 +203,12 @@ reconc_gaussian <- function(A, base_fc_mean,
 
   out <- list(
     bottom_rec_mean = as.vector(mu_b_tilde),
-    bottom_rec_covariance = Sigma_b_tilde
+    bottom_rec_cov = Sigma_b_tilde
   )
 
   if (return_upper) {
     out$upper_rec_mean <- as.vector(A %*% mu_b_tilde)
-    out$upper_rec_covariance <- A %*% Sigma_b_tilde %*% t(A)
+    out$upper_rec_cov <- A %*% Sigma_b_tilde %*% t(A)
   }
 
   return(out)

R/reconc_t.R

@@ -58,8 +58,8 @@
         )
       }
       # Seasonality test for each time series
-      is_seas = as.logical(apply(stats::ts(y_train, frequency = freq), 2, 
-                                 function(col) forecast::nsdiffs(col)))
+      is_seas = as.logical(apply(y_train, 2, 
+                                 function(col) forecast::nsdiffs(stats::ts(col,frequency = freq))))
 
     } else if (criterion == "RSS") {
       # Compute RSS for both methods, for each series
@@ -453,8 +453,7 @@ reconc_t <- function(A,
 
   out <- .core_reconc_t(
     A = A, base_fc_mean = base_fc_mean, Psi_post = Psi_post, nu_post = nu_post,
-    return_upper = return_upper, suppress_warnings = FALSE
-  )
+    return_upper = return_upper)
 
   if (return_parameters) {
     if(!is.null(posterior)) {
@@ -483,7 +482,6 @@ reconc_t <- function(A,
 #'   multivariate t-distribution.
 #' @param nu_post Degrees of freedom of the posterior multivariate t-distribution.
 #' @param return_upper Logical, whether to return the reconciled parameters for the upper variables (default is FALSE).
-#' @param suppress_warnings Logical. If TRUE, suppresses warnings about numerical issues. Default is FALSE.
 #'
 #' @return A list containing:
 #' \itemize{
@@ -501,7 +499,7 @@ reconc_t <- function(A,
 #' @keywords internal
 #' @export
 .core_reconc_t <- function(A, base_fc_mean, Psi_post, nu_post, return_upper = FALSE,
-                           return_parameters = FALSE, suppress_warnings = FALSE) {
+                           return_parameters = FALSE) {
   # Indices for Upper and Bottom
   k <- nrow(A)
   m <- ncol(A)
@@ -524,15 +522,12 @@ reconc_t <- function(A,
 
   # Invert Q using Cholesky (should be p.d.)
   Q_chol <- tryCatch(chol(Q), error = function(e) NULL)
+  # If Cholesky fails, stop
   if (is.null(Q_chol)) {
-    # Fallback to standard solve if Cholesky fails (numerical issues)
-    if (!suppress_warnings) {
-      warning("Cholesky decomposition of Q failed; using standard inversion.")
-    }
-    inv_Q <- solve(Q)
-  } else {
-    inv_Q <- chol2inv(Q_chol)
+    stop("A numerical error occured during a Cholesky decomposition. 
+         This may be due to redundant aggregated series.")
   }
+  inv_Q <- chol2inv(Q_chol)
 
   # Incoherence
   delta <- (A %*% b_hat) - u_hat