R/bounded-var.R

In jtrecenti/ea2:

###############################################################################
###############################################################################
###############################################################################

inferencia_bvar_sampling_ev <- function(bn, EV, n, delta, eps) {
  E_var <- sapply(stringr::str_match_all(as.character(EV), '(X[0-9]+) ?=='),
                  function(x) x[, 2])
  lu <- lapply(seq_along(EV), function(x) {
    range(assign_val(bn[[E_var[x]]]$cpt, EV[x])$factor)
  })
  U <- prod(sapply(lu, tail, 1))
  N_star <- 4 * log(2 / delta) * (1 + eps) / eps ^ 2
  o <- bnlearn:::schedule(bn)
  N <- M <- 0
  pa <- NULL
  X <- list()
  if (missing(n)) n <- Inf
  while (N < N_star && M < n * 100) {
    W <- 1
    X <- list()
    for (i in seq_along(o)) {
      node <- bn[[o[i]]]
      pa <- dplyr::as_data_frame(X)
      if (o[i] %in% E_var) {
        d_obs <- dplyr::filter_(node$cpt, as.character(EV[which(E_var == o[i])]))
        if (ncol(node$cpt) == 2) {
          PX <- d_obs$factor
        } else {
          nm <- intersect(names(d_obs), names(pa))
          PX <- dplyr::semi_join(d_obs, pa, nm)$factor
        }
        W <- W * PX
      } else {
        classes <- unique(node$cpt[[node$node]])
        if(ncol(node$cpt) == 2) {
          X[[o[i]]] <- classes[LaplacesDemonCpp::rcat(1, node$cpt$factor)]
        } else {
          nm <- intersect(names(node$cpt), names(pa))
          probs <- dplyr::semi_join(node$cpt, pa, nm)$factor
          X[[o[i]]] <- classes[LaplacesDemonCpp::rcat(1, probs)]
        }
      }
    }
    N <- N + W / U
    M <- M + 1
  }
  U * N / M
}

inferencia_bvar_sampling <- function(bn, Q, E, n, delta = 0.9, eps = 0.1) {
  if (missing(Q) && missing(E)) return(NULL)
  if (!missing(Q) && !missing(E)) {
    EV <- c(Q, E)
    RAZ <- E
    p1 <- inferencia_bvar_sampling_ev(bn, EV, n, delta, eps)
    p2 <- inferencia_bvar_sampling_ev(bn, RAZ, n, delta, eps)
    return(p1 / p2)
  } else if (missing(Q) | !missing(E)) {
    EV <- E
  } else {
    EV <- Q
  }
  return(inferencia_bvar_sampling_ev(bn, EV, n, delta, eps))
}