R/pmvt.R

Defines functions pmvt

Documented in pmvt

#' Compute multivariate Student-t (MVT) probabilities that have spatial covariance
#' matrices using Vecchia approximation
#' @importFrom truncnorm etruncnorm
#' @importFrom nleqslv nleqslv
#' @importFrom TruncatedNormal cholperm
#' @importFrom utils getFromNamespace
#'
#' @param lower lower bound vector for TMVT
#' @param upper upper bound vector for TMVT
#' @param delta MVT shifting parameter
#' @param df degrees of freedom
#' @param locs location (feature) matrix n X d
#' @param covName covariance function name from the `GpGp` package
#' @param covParms parameters for `covName`
#' @param m Vecchia conditioning set size
#' @param sigma dense covariance matrix, not needed when `locs` is not null
#' @param reorder whether to reorder integration variables. `0` for no,
#' `1` for FIC-based univariate ordering, `2` for Vecchia-based univariate
#' ordering, and `3` for univariate reordering, which appeared faster than `2`
#' @param NLevel1 first level Monte Carlo sample size
#' @param NLevel2 second level Monte Carlo sample size
#' @param verbose verbose or not
#' @param retlog TRUE or FALSE for whether to return loglk or not
#' @param ... could be
#' m_ord for conditioning set size for reordering
#' @return estimated MVT probability and estimation error
#'
#' @export
pmvt <- function(lower, upper, delta, df, locs = NULL, covName = "matern15_isotropic",
                 covParms = c(1.0, 0.1, 0.0), m = 30, sigma = NULL, reorder = 0,
                 NLevel1 = 12, NLevel2 = 1e4, verbose = FALSE, retlog = FALSE, ...) {
  # standardize the input MVN prob -----------------------------
  lower <- lower - delta
  upper <- upper - delta
  if (is.null(sigma)) {
    n <- nrow(locs)
    use_sigma <- FALSE
    margin_sd <- sqrt(covParms[1])
    upper <- upper / margin_sd
    lower <- lower / margin_sd
    covParms[1] <- 1
  } else {
    n <- nrow(sigma)
    use_sigma <- TRUE
    margin_sd <- sqrt(diag(sigma))
    upper <- upper / margin_sd
    lower <- lower / margin_sd
    sigma <- t(t(sigma / margin_sd) / margin_sd)
  }
  if (any(upper < lower)) {
    stop("Invalid MVN probability. Truncated marginal
         probabilities have negative value(s)\n")
  }
  lower_upper <- matrix(0, n, 2)
  lower_upper[, 1] <- lower
  lower_upper[, 2] <- upper
  lower <- lower_upper[, 1]
  upper <- lower_upper[, 2]
  # reorder --------------------------------
  if (is.null(list(...)[["m_ord"]])) {
    m_ord <- m
  } else {
    m_ord <- list(...)[["m_ord"]]
  }
  if (reorder == 1) {
    if (use_sigma) {
      ord <- FIC_reorder_univar(lower, upper, m_ord, covMat = sigma)
    } else {
      ord <- FIC_reorder_univar(
        lower, upper, m_ord, locs, covName,
        covParms
      )
    }
    lower <- lower[ord]
    upper <- upper[ord]
    if (use_sigma) {
      sigma <- sigma[ord, ord, drop = FALSE]
    } else {
      locs <- locs[ord, , drop = FALSE]
    }
  } else if (reorder == 2) {
    if (use_sigma) {
      ord <- Vecc_reorder(lower, upper, m_ord, covMat = sigma)$order
    } else {
      ord <- Vecc_reorder(
        lower, upper, m_ord, locs, covName, covParms
      )$order
    }
    lower <- lower[ord]
    upper <- upper[ord]
    if (use_sigma) {
      sigma <- sigma[ord, ord, drop = FALSE]
    } else {
      locs <- locs[ord, , drop = FALSE]
    }
  } else if (reorder == 3) {
    if (!use_sigma) {
      cov_func_GpGp <- utils::getFromNamespace(covName, "GpGp")
      sigma <- cov_func_GpGp(covParms, locs)
    }
    ord <- univar_order(lower, upper, sigma)
    lower <- lower[ord]
    upper <- upper[ord]
    if (use_sigma) {
      sigma <- sigma[ord, ord, drop = FALSE]
    } else {
      locs <- locs[ord, , drop = FALSE]
    }
  }
  # find nearest neighbors for Vecchia --------------------------------
  if (use_sigma) {
    NN <- find_nn_corr(sigma, m)
  } else {
    NN <- GpGp::find_ordered_nn(locs, m)
  }
  # find Vecchia approx object -----------------------------------
  if (use_sigma) {
    vecc_obj <- vecc_cond_mean_var_sp(NN, covMat = sigma)
    idx <- which(vecc_obj$cond_var < 0.01)
    if (length(idx) > 0) {
      diag(sigma)[idx] <- diag(sigma)[idx] + 0.01
      vecc_obj <- vecc_cond_mean_var_sp(NN, covMat = sigma)
    }
  } else {
    vecc_obj <- vecc_cond_mean_var_sp(NN,
      locs = locs, covName = covName,
      covParms = covParms
    )
    if (any(vecc_obj$cond_var < 0.01)) {
      covParms[length(covParms)] <- 0.01 # nugget
      vecc_obj <- vecc_cond_mean_var_sp(NN,
        locs = locs, covName = covName,
        covParms = covParms
      )
    }
  }
  # find tilting parameter beta -----------------------------------
  trunc_expect <- truncnorm::etruncnorm(lower, upper)
  x0 <- c(trunc_expect, rep(0, n))
  x0[2 * n] <- sqrt(df)
  x0[n] <- log(x0[2 * n])
  solv_idea_5_sp <- nleqslv::nleqslv(
    x = x0, fn = gradpsiT_idea5, veccCondMeanVarObj = vecc_obj,
    a = lower, b = upper, nu = df, global = "pwldog", method = "Broyden"
  )
  soln <- solv_idea_5_sp$x
  exitflag <- solv_idea_5_sp$termcd
  if(!(exitflag %in% c(1,2)) || !all.equal(solv_idea_5_sp$fvec, rep(0, length(x0)))){
    warning('Did not find a solution to the nonlinear system in `pmvt`!')
  }
  soln[n] <- exp(soln[n])
  x <- soln[1:n]
  beta <- soln[(n + 1):(2 * n)]
  eta <- beta[n]
  # compute MVT probs and est error ---------------------------------
  const <- log(2*pi) / 2 - lgamma(df / 2) - (df / 2 - 1) * log(2) +
    TruncatedNormal::lnNpr(-eta, Inf) + 0.5 * eta^2;
  exp_psi <- sample_psiT_idea5_cpp(vecc_obj, lower, upper, df,
    beta = beta, N_level1 = NLevel1,
    N_level2 = NLevel2
  )
  if (retlog) {
    exponent <- min(exp_psi[[2]])
    log_est_prob <- exponent +
      log(mean(exp_psi[[1]] * exp(exp_psi[[2]] - exponent))) + const
    return(log_est_prob)
  } else {
    exp_psi <- exp_psi[[1]] * exp(exp_psi[[2]])
    est_prob <- mean(exp_psi) * exp(const)
    est_prob_err <- stats::sd(exp_psi) / sqrt(NLevel1) * exp(const)
    attr(est_prob, "error") <- est_prob_err
    return(est_prob)
  }
}


# TEST -------------------------------------------------------

# library(GpGp)
# library(TruncatedNormal)
# library(mvtnorm)
# library(VeccTMVN)
# 
# 
# ## example MVN probabilities --------------------------------
# set.seed(123)
# n1 <- 10
# n2 <- 10
# n <- n1 * n2
# locs <- as.matrix(expand.grid((1:n1) / n1, (1:n2) / n2))
# covparms <- c(2, 0.3, 0)
# nu <- 10
# cov_mat <- matern15_isotropic(covparms, locs)
# a_list <- list(rep(-Inf, n), rep(-1, n), -runif(n) * 2 - 4)
# b_list <- list(rep(-2, n), rep(1, n), -runif(n) * 2)
# 
# ## Compute MVN probs --------------------------------
# N_level1 <- 12 # Level 1 MC size
# N_level2 <- 1e4 # Level 2 MC size
# m <- 30 # num of nearest neighbors
# for (i in 1:length(a_list)) {
#   a <- a_list[[i]]
#   b <- b_list[[i]]
#   
#   VeccTMVN::pmvt(a, b, 0, nu, locs,
#                  covName = "matern15_isotropic",
#                  covParms = covparms, m = m, verbose = FALSE
#   )
#   
#   ### Compute MVN prob with idea V -----------------------
#   time_Vecc <- system.time(est_Vecc <- VeccTMVN::pmvt(a, b, 0, nu, locs,
#                                                       covName = "matern15_isotropic",
#                                                       covParms = covparms, m = m, verbose = FALSE
#   ))[[3]]
#   ### Compute MVN prob with other methods -----------------------
#   time_TN <- system.time(est_TN <- TruncatedNormal::pmvt(
#     rep(0, n), cov_mat,
#     lb = a, ub = b, df = nu
#   ))[[3]]
#   time_TLR <- system.time(
#     est_TLR <- tlrmvnmvt::pmvt(a, b, df = nu, sigma = cov_mat)
#   )[[3]]
#   time_Genz <- system.time(
#     est_Genz <- mvtnorm::pmvt(a, b, df = nu, sigma = cov_mat)
#   )[[3]]
#   cat(
#     "est_Vecc", est_Vecc, "err_Vecc", attributes(est_Vecc)$error, "time_Vecc",
#     time_Vecc, "\n"
#   )
#   cat(
#     "est_TN", est_TN, "err_TN", attributes(est_TN)$relerr * est_TN, "time_TN",
#     time_TN, "\n"
#   )
#   cat(
#     "est_TLR", est_TLR, "err_TLR", attributes(est_TLR)$error, "time_TLR",
#     time_TLR, "\n"
#   )
#   cat(
#     "est_Genz", est_Genz, "err_Genz", attributes(est_Genz)$error, "time_Genz",
#     time_Genz, "\n"
#   )
# }

Try the VeccTMVN package in your browser

Any scripts or data that you put into this service are public.

VeccTMVN documentation built on Aug. 21, 2025, 6:01 p.m.