restructure bootstrap, intervals for coefs and probs

helske · Oct 11, 2024 · 5f58906 · 5f58906
1 parent 24f01ab
commit 5f58906
Show file tree

Hide file tree

Showing 14 changed files with 426 additions and 490 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -101,14 +101,6 @@ get_log_pi <- function(gamma, X) {
     .Call(`_seqHMM_get_log_pi`, gamma, X)
 }
 
-get_pi_all <- function(gamma, X) {
-    .Call(`_seqHMM_get_pi_all`, gamma, X)
-}
-
-get_pi_boot <- function(gamma, X) {
-    .Call(`_seqHMM_get_pi_boot`, gamma, X)
-}
-
 get_A <- function(gamma, X, tv) {
     .Call(`_seqHMM_get_A`, gamma, X, tv)
 }
@@ -117,10 +109,6 @@ get_log_A <- function(gamma, X, tv) {
     .Call(`_seqHMM_get_log_A`, gamma, X, tv)
 }
 
-get_A_all <- function(gamma, X, tv) {
-    .Call(`_seqHMM_get_A_all`, gamma, X, tv)
-}
-
 get_B <- function(gamma, X, add_missing, tv) {
     .Call(`_seqHMM_get_B`, gamma, X, add_missing, tv)
 }
@@ -129,10 +117,34 @@ get_log_B <- function(gamma, X, add_missing, tv) {
     .Call(`_seqHMM_get_log_B`, gamma, X, add_missing, tv)
 }
 
+get_pi_all <- function(gamma, X) {
+    .Call(`_seqHMM_get_pi_all`, gamma, X)
+}
+
+get_A_all <- function(gamma, X, tv) {
+    .Call(`_seqHMM_get_A_all`, gamma, X, tv)
+}
+
 get_B_all <- function(gamma, X, add_missing, tv) {
     .Call(`_seqHMM_get_B_all`, gamma, X, add_missing, tv)
 }
 
+get_pi_boot <- function(gamma, X) {
+    .Call(`_seqHMM_get_pi_boot`, gamma, X)
+}
+
+get_A_boot <- function(gamma, X, tv) {
+    .Call(`_seqHMM_get_A_boot`, gamma, X, tv)
+}
+
+get_B_boot <- function(gamma, X, tv) {
+    .Call(`_seqHMM_get_B_boot`, gamma, X, tv)
+}
+
+get_omega_boot <- function(gamma, X) {
+    .Call(`_seqHMM_get_omega_boot`, gamma, X)
+}
+
 logLikHMM <- function(transition, emission, init, obs, threads) {
     .Call(`_seqHMM_logLikHMM`, transition, emission, init, obs, threads)
 }

diff --git a/R/bootstrap.R b/R/bootstrap.R
@@ -78,13 +78,18 @@ bootstrap_coefs.nhmm <- function(model, B = 1000,
     penalty <- model$estimation_results$penalty
   }
   gammas_mle <- model$gammas
-  coefs <- matrix(NA, length(unlist(gammas_mle)), B)
+  gamma_pi <- replicate(B, gammas_mle$pi, simplify = FALSE)
+  gamma_A <- replicate(B, gammas_mle$A, simplify = FALSE)
+  gamma_B <- replicate(B, gammas_mle$B, simplify = FALSE)
   if (verbose) pb <- utils::txtProgressBar(min = 0, max = 100, style = 3)
   if (method == "nonparametric") {
     for (i in seq_len(B)) {
       mod <- bootstrap_model(model)
       fit <- fit_nhmm(mod, init, 0, 0, 1, penalty, ...)
-      coefs[, i] <- unlist(permute_states(fit$gammas, gammas_mle))
+      fit$gammas <- permute_states(fit$gammas, gammas_mle)
+      gamma_pi[[i]] <- fit$gammas$pi
+      gamma_A[[i]] <- fit$gammas$A
+      gamma_B[[i]] <- fit$gammas$B
       if (verbose) utils::setTxtProgressBar(pb, 100 * i/B)
     }
   } else {
@@ -103,12 +108,15 @@ bootstrap_coefs.nhmm <- function(model, B = 1000,
         N, T_, M, S, formula_pi, formula_A, formula_B,
         data = d, time, id, init)$model
       fit <- fit_nhmm(mod, init, 0, 0, 1, penalty, ...)
-      coefs[, i] <- unlist(permute_states(fit$gammas, gammas_mle))
+      fit$gammas <- permute_states(fit$gammas, gammas_mle)
+      gamma_pi[[i]] <- fit$gammas$pi
+      gamma_A[[i]] <- fit$gammas$A
+      gamma_B[[i]] <- fit$gammas$B
       if (verbose) utils::setTxtProgressBar(pb, 100 * i/B)
     }
   }
   if (verbose) close(pb)
-  model$boot <- coefs
+  model$boot <- list(gamma_pi = gamma_pi, gamma_A = gamma_A, gamma_B = gamma_B)
   model
 }
 #' @rdname bootstrap
@@ -126,7 +134,10 @@ bootstrap_coefs.mnhmm <- function(model, B = 1000,
     penalty <- model$estimation_results$penalty
   }
   gammas_mle <- model$gammas
-  coefs <- matrix(NA, length(unlist(gammas_mle)), B)
+  gamma_pi <- replicate(B, gammas_mle$pi, simplify = FALSE)
+  gamma_A <- replicate(B, gammas_mle$A, simplify = FALSE)
+  gamma_B <- replicate(B, gammas_mle$B, simplify = FALSE)
+  gamma_omega <- replicate(B, gammas_mle$omega, simplify = FALSE)
   D <- model$n_clusters
   if (verbose) pb <- utils::txtProgressBar(min = 0, max = 100, style = 3)
   if (method == "nonparametric") {
@@ -142,7 +153,10 @@ bootstrap_coefs.mnhmm <- function(model, B = 1000,
         fit$gammas$A[[j]] <- out$A
         fit$gammas$B[[j]] <- out$B
       }
-      coefs[, i] <- unlist(fit$gammas)
+      gamma_pi[[i]] <- fit$gammas$pi
+      gamma_A[[i]] <- fit$gammas$A
+      gamma_B[[i]] <- fit$gammas$B
+      gamma_omega[[i]] <- fit$gammas$omega
       if (verbose) utils::setTxtProgressBar(pb, 100 * i/B)
     }
   } else {
@@ -171,11 +185,17 @@ bootstrap_coefs.mnhmm <- function(model, B = 1000,
         fit$gammas$A[[j]] <- out$A
         fit$gammas$B[[j]] <- out$B
       }
-      coefs[, i] <- unlist(fit$gammas)
+      gamma_pi[[i]] <- fit$gammas$pi
+      gamma_A[[i]] <- fit$gammas$A
+      gamma_B[[i]] <- fit$gammas$B
+      gamma_omega[[i]] <- fit$gammas$omega
       if (verbose) utils::setTxtProgressBar(pb, 100 * i/B)
     }
   }
   if (verbose) close(pb)
-  model$boot <- coefs
+  model$boot <- list(
+    gamma_pi = gamma_pi, gamma_A = gamma_A, gamma_B = gamma_B,
+    gamma_omega = gamma_omega
+  )
   model
 }
diff --git a/R/coef.R b/R/coef.R
@@ -47,21 +47,24 @@ coef.nhmm <- function(object, probs, ...) {
       "Argument {.arg probs} must be a {.cls numeric} vector with values
       between 0 and 1."
     )
-    p_i <- length(gamma_pi)
-    p_s <- length(gamma_A)
-    p_o <- length(gamma_B)
+    p_i <- length(unlist(object$gammas$pi))
+    p_s <- length(unlist(object$gammas$A))
+    p_o <- length(unlist(object$gammas$B))
     stopifnot_(
       !is.null(object$boot),
       paste0(
         "Model does not contain bootstrap samples of coefficients. ",
         "Run {.fn bootstrap_coefs} first."
       )
     )
-    qs <- seqHMM:::fast_quantiles(object$boot, probs)
+    B <- length(object$boot$gamma_pi)
+    q_pi <- seqHMM:::fast_quantiles(matrix(unlist(object$boot$gamma_pi), ncol = B), probs)
+    q_A <- seqHMM:::fast_quantiles(matrix(unlist(object$boot$gamma_A), ncol = B), probs)
+    q_B <- seqHMM:::fast_quantiles(matrix(unlist(object$boot$gamma_B), ncol = B), probs)
     for(i in seq_along(probs)) {
-      gamma_pi[paste0("q", 100 * probs[i])] <- qs[seq_len(p_i), i]
-      gamma_A[paste0("q", 100 * probs[i])] <- qs[p_i + seq_len(p_s), i]
-      gamma_B[paste0("q", 100 * probs[i])] <- qs[p_i + p_s + seq_len(p_o), i]
+      gamma_pi[paste0("q", 100 * probs[i])] <- q_pi[, i]
+      gamma_A[paste0("q", 100 * probs[i])] <- q_A[, i]
+      gamma_B[paste0("q", 100 * probs[i])] <- q_B[, i]
     }
   }
   list(
@@ -80,39 +83,39 @@ coef.mnhmm <- function(object, probs, ...) {
   K_i <- length(object$coef_names_initial)
   object$state_names <- unname(object$state_names)
   gamma_pi <- data.frame(
+    cluster = rep(object$cluster_names, each = S * K_i),
     state = unlist(object$state_names),
     parameter = rep(object$coef_names_initial, each = S),
-    estimate = unlist(object$gammas$pi),
-    cluster = rep(object$cluster_names, each = S * K_i)
+    estimate = unlist(object$gammas$pi)
   )
   K_s <- length(object$coef_names_transition)
   gamma_A <- data.frame(
+    cluster = rep(object$cluster_names, each = S * S * K_s),
     state_from = unlist(object$state_names),
     state_to = rep(unlist(object$state_names), each = S),
     parameter = rep(object$coef_names_transition, each = S * S),
-    estimate = unlist(object$gammas$A),
-    cluster = rep(object$cluster_names, each = S * S * K_s)
+    estimate = unlist(object$gammas$A)
   )
   K_o <- length(object$coef_names_emission)
   if (object$n_channels == 1) {
     gamma_B <- data.frame(
+      cluster =  rep(object$cluster_names, each = S * M * K_o),
       state = unlist(object$state_names),
       observations = rep(object$symbol_names, each = S),
       parameter = rep(object$coef_names_emission, each = S * M),
-      estimate = unlist(object$gammas$B),
-      cluster =  rep(object$cluster_names, each = S * M * K_o)
+      estimate = unlist(object$gammas$B)
     )
   } else {
     gamma_B <- data.frame(
+      cluster =  unlist(lapply(seq_len(object$n_channels), function(i) {
+        rep(object$cluster_names, each = S * M[i] * K_o)
+      })),
       state = unlist(object$state_names),
       observations = rep(unlist(object$symbol_names), each = S),
       parameter = unlist(lapply(seq_len(object$n_channels), function(i) {
         rep(object$coef_names_emission, each = S * M[i])
       })),
-      estimate = unlist(object$gammas$B),
-      cluster =  unlist(lapply(seq_len(object$n_channels), function(i) {
-        rep(object$cluster_names, each = S * M[i] * K_o)
-      }))
+      estimate = unlist(object$gammas$B)
     )
   }
   gamma_omega <- data.frame(
@@ -128,10 +131,10 @@ coef.mnhmm <- function(object, probs, ...) {
       "Argument {.arg probs} must be a {.cls numeric} vector with values
       between 0 and 1."
     )
-    p_i <- length(gamma_pi)
-    p_s <- length(gamma_A)
-    p_o <- length(gamma_B)
-    p_d <- length(gamma_omega)
+    p_i <- length(unlist(object$gammas$pi))
+    p_s <- length(unlist(object$gammas$A))
+    p_o <- length(unlist(object$gammas$B))
+    p_d <- length(object$gammas$omega)
 
     stopifnot_(
       !is.null(object$boot),
@@ -140,13 +143,16 @@ coef.mnhmm <- function(object, probs, ...) {
         "Run {.fn bootstrap_coefs} first."
       )
     )
-    qs <- seqHMM:::fast_quantiles(object$boot, probs)
+    B <- length(object$boot$gamma_pi)
+    q_pi <- seqHMM:::fast_quantiles(matrix(unlist(object$boot$gamma_pi), ncol = B), probs)
+    q_A <- seqHMM:::fast_quantiles(matrix(unlist(object$boot$gamma_A), ncol = B), probs)
+    q_B <- seqHMM:::fast_quantiles(matrix(unlist(object$boot$gamma_B), ncol = B), probs)
+    q_omega <- seqHMM:::fast_quantiles(matrix(unlist(object$boot$gamma_omega), ncol = B), probs)
     for(i in seq_along(probs)) {
-      gamma_pi[paste0("q", 100 * probs[i])] <- qs[seq_len(p_i), i]
-      gamma_A[paste0("q", 100 * probs[i])] <- qs[p_i + seq_len(p_s), i]
-      gamma_B[paste0("q", 100 * probs[i])] <- qs[p_i + p_s + seq_len(p_o), i]
-      gamma_omega[paste0("q", 100 * probs[i])] <- 
-        qs[p_i + p_s + p_o + seq_len(p_d), i]
+      gamma_pi[paste0("q", 100 * probs[i])] <- q_pi[, i]
+      gamma_A[paste0("q", 100 * probs[i])] <- q_A[, i]
+      gamma_B[paste0("q", 100 * probs[i])] <- q_B[, i]
+      gamma_omega[paste0("q", 100 * probs[i])] <- q_omega[, i]
     }
   }
   list(

diff --git a/R/create_initial_values.R b/R/create_initial_values.R
@@ -123,16 +123,16 @@ create_eta_B_inits <- function(x, S, M, K, init_sd = 0, D = 1) {
       if (is.null(x)) {
         create_eta_multichannel_B_nhmm(
           numeric(sum((M - 1) * K * S)), S, M, K, init_sd
-          )
+        )
       } else {
         stopifnot_(
-          length(x) == sum((M - 1) * K * S),
+          length(unlist(x)) == sum((M - 1) * K * S),
           paste0(
             "Number of initial values for {.val eta_B} is not equal to ",
             "sum((M - 1) * K * S) = {sum((M - 1) * K * S)}."
           )
         )
-        create_eta_multichannel_B_nhmm(x, S, M, K, init_sd)
+        create_eta_multichannel_B_nhmm(unlist(x), S, M, K, init_sd)
       }
     }
   } else {

diff --git a/R/fit_nhmm.R b/R/fit_nhmm.R
@@ -77,9 +77,6 @@ fit_nhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
     Qm <- t(create_Q(M))
     if (need_grad) {
       objectivef <- function(pars) {
-        # if (any(!is.finite(exp(pars)))) {
-        #   return(list(objective = Inf, gradient = rep(-Inf, length(pars))))
-        # }
         eta_pi <- create_eta_pi_nhmm(pars[seq_len(n_i)], S, K_i)
         eta_A <- create_eta_A_nhmm(pars[n_i + seq_len(n_s)], S, K_s)
         eta_B <- create_eta_B_nhmm(
@@ -94,9 +91,6 @@ fit_nhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
       }
     } else {
       objectivef <- function(pars) {
-        # if (any(!is.finite(exp(pars)))) {
-        #   return(Inf)
-        # } 
         eta_pi <- create_eta_pi_nhmm(pars[seq_len(n_i)], S, K_i)
         eta_A <- create_eta_A_nhmm(pars[n_i + seq_len(n_s)], S, K_s)
         eta_B <- create_eta_B_nhmm(
@@ -114,9 +108,6 @@ fit_nhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
     Qm <- lapply(M, function(m) t(create_Q(m)))
     if (need_grad) {
       objectivef <- function(pars) {
-        # if (any(!is.finite(exp(pars)))) {
-        #   return(list(objective = Inf, gradient = rep(-Inf, length(pars))))
-        # }
         eta_pi <- create_eta_pi_nhmm(pars[seq_len(n_i)], S, K_i)
         eta_A <- create_eta_A_nhmm(pars[n_i + seq_len(n_s)], S, K_s)
         eta_B <- create_eta_multichannel_B_nhmm(
@@ -131,9 +122,6 @@ fit_nhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
       }
     } else {
       objectivef <- function(pars) {
-        # if (any(!is.finite(exp(pars)))) {
-        #   return(Inf)
-        # }
         eta_pi <- create_eta_pi_nhmm(pars[seq_len(n_i)], S, K_i)
         eta_A <- create_eta_A_nhmm(pars[n_i + seq_len(n_s)], S, K_s)
         eta_B <- create_eta_multichannel_B_nhmm(