get_probs, new methods, fix c++ indices, test gradients

helske · Sep 23, 2024 · bc57af3 · bc57af3
1 parent 709658f
commit bc57af3
Show file tree

Hide file tree

Showing 40 changed files with 1,712 additions and 719 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -4,6 +4,8 @@ S3method("cluster_names<-",mhmm)
 S3method("cluster_names<-",mnhmm)
 S3method("state_names<-",hmm)
 S3method("state_names<-",mhmm)
+S3method("state_names<-",mnhmm)
+S3method("state_names<-",nhmm)
 S3method(ame,mnhmm)
 S3method(ame,nhmm)
 S3method(cluster_names,mhmm)
@@ -14,8 +16,22 @@ S3method(forward_backward,hmm)
 S3method(forward_backward,mhmm)
 S3method(forward_backward,mnhmm)
 S3method(forward_backward,nhmm)
+S3method(get_cluster_probs,mhmm)
+S3method(get_cluster_probs,mnhmm)
+S3method(get_emission_probs,hmm)
+S3method(get_emission_probs,mhmm)
+S3method(get_emission_probs,mnhmm)
+S3method(get_emission_probs,nhmm)
+S3method(get_initial_probs,hmm)
+S3method(get_initial_probs,mhmm)
+S3method(get_initial_probs,mnhmm)
+S3method(get_initial_probs,nhmm)
 S3method(get_probs,mnhmm)
 S3method(get_probs,nhmm)
+S3method(get_transition_probs,hmm)
+S3method(get_transition_probs,mhmm)
+S3method(get_transition_probs,mnhmm)
+S3method(get_transition_probs,nhmm)
 S3method(hidden_paths,hmm)
 S3method(hidden_paths,mhmm)
 S3method(hidden_paths,mnhmm)
@@ -41,6 +57,8 @@ S3method(print,summary.mnhmm)
 S3method(print,summary.nhmm)
 S3method(state_names,hmm)
 S3method(state_names,mhmm)
+S3method(state_names,mnhmm)
+S3method(state_names,nhmm)
 S3method(summary,mhmm)
 S3method(summary,mnhmm)
 S3method(summary,nhmm)
@@ -61,7 +79,11 @@ export(estimate_mnhmm)
 export(estimate_nhmm)
 export(fit_model)
 export(forward_backward)
+export(get_cluster_probs)
+export(get_emission_probs)
+export(get_initial_probs)
 export(get_probs)
+export(get_transition_probs)
 export(gridplot)
 export(hidden_paths)
 export(mc_to_sc)

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -61,12 +61,12 @@ get_pi <- function(gamma_raw, X, logspace) {
     .Call(`_seqHMM_get_pi`, gamma_raw, X, logspace)
 }
 
-get_A <- function(gamma_raw, X, logspace) {
-    .Call(`_seqHMM_get_A`, gamma_raw, X, logspace)
+get_A <- function(gamma_raw, X, logspace, tv) {
+    .Call(`_seqHMM_get_A`, gamma_raw, X, logspace, tv)
 }
 
-get_B <- function(gamma_raw, X, M, logspace, add_missing) {
-    .Call(`_seqHMM_get_B`, gamma_raw, X, M, logspace, add_missing)
+get_B <- function(gamma_raw, X, M, logspace, add_missing, tv) {
+    .Call(`_seqHMM_get_B`, gamma_raw, X, M, logspace, add_missing, tv)
 }
 
 logLikHMM <- function(transition, emission, init, obs, threads) {
@@ -109,20 +109,20 @@ log_objective <- function(transition, emission, init, obs, ANZ, BNZ, INZ, nSymbo
     .Call(`_seqHMM_log_objective`, transition, emission, init, obs, ANZ, BNZ, INZ, nSymbols, threads)
 }
 
-log_objective_nhmm_singlechannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs) {
-    .Call(`_seqHMM_log_objective_nhmm_singlechannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs)
+log_objective_nhmm_singlechannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, iv_pi, iv_A, iv_B, tv_A, tv_B) {
+    .Call(`_seqHMM_log_objective_nhmm_singlechannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, iv_pi, iv_A, iv_B, tv_A, tv_B)
 }
 
-log_objective_nhmm_multichannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, M) {
-    .Call(`_seqHMM_log_objective_nhmm_multichannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, M)
+log_objective_nhmm_multichannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, M, iv_pi, iv_A, iv_B, tv_A, tv_B) {
+    .Call(`_seqHMM_log_objective_nhmm_multichannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, M, iv_pi, iv_A, iv_B, tv_A, tv_B)
 }
 
-log_objective_mnhmm_singlechannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs) {
-    .Call(`_seqHMM_log_objective_mnhmm_singlechannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs)
+log_objective_mnhmm_singlechannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs, iv_pi, iv_A, iv_B, tv_A, tv_B, iv_omega) {
+    .Call(`_seqHMM_log_objective_mnhmm_singlechannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs, iv_pi, iv_A, iv_B, tv_A, tv_B, iv_omega)
 }
 
-log_objective_mnhmm_multichannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs, M) {
-    .Call(`_seqHMM_log_objective_mnhmm_multichannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs, M)
+log_objective_mnhmm_multichannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs, M, iv_pi, iv_A, iv_B, tv_A, tv_B, iv_omega) {
+    .Call(`_seqHMM_log_objective_mnhmm_multichannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs, M, iv_pi, iv_A, iv_B, tv_A, tv_B, iv_omega)
 }
 
 log_objectivex <- function(transition, emission, init, obs, ANZ, BNZ, INZ, nSymbols, coef, X, numberOfStates, threads) {

diff --git a/R/ame.R b/R/ame.R
@@ -72,9 +72,9 @@ ame.nhmm <- function(
   }
   # use same RNG seed so that the same samples of coefficients are drawn
   newdata[[variable]] <- values[1]
-  pred <- predict(model, newdata, dontchange_colnames = TRUE)
+  pred <- get_probs(model, newdata)
   newdata[[variable]] <- values[2]
-  pred2 <- predict(model, newdata, dontchange_colnames = TRUE)
+  pred2 <- get_probs(model, newdata)
   pars <- c("initial_probs", "transition_probs", "emission_probs")
   for (i in pars) {
     pred[[i]]$estimate <- pred[[i]]$estimate - pred2[[i]]$estimate
@@ -170,10 +170,10 @@ ame.mnhmm <- function(
   }
   # use same RNG seed so that the same samples of coefficients are drawn
   newdata[[variable]] <- values[1]
-  pred <- predict(model, newdata, dontchange_colnames = TRUE)
+  pred <- get_probs(model, newdata)
 
   newdata[[variable]] <- values[2]
-  pred2 <- predict(model, newdata, dontchange_colnames = TRUE)
+  pred2 <- get_probs(model, newdata)
   pars <- c("initial_probs", "transition_probs", "emission_probs",
             "cluster_probs")
   for (i in pars) {
@@ -186,18 +186,18 @@ ame.mnhmm <- function(
     "states" = c("cluster", "time", channel, "observation"),
     "sequences" = c("cluster", "time", "state", channel, "observation")
   )
-
+  
   out <- list()
   out$initial_probs <- pred$initial_probs |> 
-      dplyr::group_by(cluster, state) |>
-      dplyr::summarise(estimate = mean(estimate)) |> 
-      dplyr::ungroup()
+    dplyr::group_by(cluster, state) |>
+    dplyr::summarise(estimate = mean(estimate)) |> 
+    dplyr::ungroup()
 
   out$transition_probs <- pred$transition_probs |> 
-      dplyr::group_by(cluster, time, state_from, state_to) |>
-      dplyr::summarise(estimate = mean(estimate)) |> 
-      dplyr::ungroup() |> 
-      dplyr::rename(!!time := time)
+    dplyr::group_by(cluster, time, state_from, state_to) |>
+    dplyr::summarise(estimate = mean(estimate)) |> 
+    dplyr::ungroup() |> 
+    dplyr::rename(!!time := time)
 
   out$emission_probs <- pred$emission_probs |> 
     dplyr::group_by(dplyr::pick(group_by_B)) |>

diff --git a/R/build_mnhmm.R b/R/build_mnhmm.R
@@ -27,8 +27,19 @@ build_mnhmm <- function(
     out$model,
     class = "mnhmm",
     nobs = attr(out$observations, "nobs"),
-    df = out$extras$n_pars,
+    df = out$extras$np_omega + 
+      n_clusters * (out$extras$np_pi + out$extras$np_A + out$extras$np_B),
     type = paste0(out$extras$multichannel, "mnhmm"),
-    intercept_only = out$extras$intercept_only
+    intercept_only = out$extras$intercept_only,
+    iv_pi = out$extras$iv_pi,
+    iv_A =  out$extras$iv_A,
+    iv_B = out$extras$iv_B,
+    iv_omega =  out$extras$iv_omega,
+    tv_A =  out$extras$tv_A,
+    tv_B =  out$extras$tv_B,
+    np_pi = n_clusters * out$extras$np_pi,
+    np_A = n_clusters * out$extras$np_A,
+    np_B = n_clusters * out$extras$np_B,
+    np_omega = out$extras$np_omega
   )
 }
diff --git a/R/build_nhmm.R b/R/build_nhmm.R
@@ -14,8 +14,16 @@ build_nhmm <- function(
     out$model,
     class = "nhmm",
     nobs = attr(out$observations, "nobs"),
-    df = out$extras$n_pars,
+    df = out$extras$np_pi + out$extras$np_A + out$extras$np_B,
     type = paste0(out$extras$multichannel, "nhmm"),
-    intercept_only = out$extras$intercept_only
+    intercept_only = out$extras$intercept_only,
+    iv_pi = out$extras$iv_pi,
+    iv_A =  out$extras$iv_A,
+    iv_B = out$extras$iv_B,
+    tv_A =  out$extras$tv_A,
+    tv_B =  out$extras$tv_B,
+    np_pi = out$extras$np_pi,
+    np_A = out$extras$np_A,
+    np_B = out$extras$np_B
   )
 }
diff --git a/R/coef.R b/R/coef.R
@@ -43,47 +43,47 @@ coef.nhmm <- function(object, probs = c(0.025, 0.975), ...) {
     )
   }
 
-  stopifnot_(
-    checkmate::test_numeric(
-      x = probs, lower = 0, upper = 1, any.missing = FALSE, min.len = 1L
-    ),
-    "Argument {.arg probs} must be a {.cls numeric} vector with values
-      between 0 and 1."
-  )
-  p_i <- length(gamma_pi_raw)
-  p_s <- length(gamma_A_raw)
-  p_o <- length(gamma_B_raw)
-  sds <- try(
-    diag(solve(-object$estimation_results$hessian)), 
-    silent = TRUE
-  )
-  if (inherits(sds, "try-error")) {
-    warning_(
-      paste0(
-        "Standard errors could not be computed due to singular Hessian. ",
-        "Confidence intervals will not be provided."
-      )
-    )
-    sds <- rep(NA, p_i + p_s + p_o)
-  } else {
-    if (any(sds < 0)) {
-      warning_(
-        paste0(
-          "Standard errors could not be computed due to negative variances. ",
-          "Confidence intervals will not be provided."
-        )
-      )
-      sds <- rep(NA, p_i + p_s + p_o)
-    } else {
-      sds <- sqrt(sds)
-    }
-  }
-  for(i in seq_along(probs)) {
-    q <- qnorm(probs[i])
-    gamma_pi[paste0("q", 100 * probs[i])] <- gamma_pi_raw + q * sds[seq_len(p_i)]
-    gamma_A[paste0("q", 100 * probs[i])] <- gamma_A_raw + q * sds[p_i + seq_len(p_s)]
-    gamma_B[paste0("q", 100 * probs[i])] <- gamma_B_raw + q * sds[p_i + p_s + seq_len(p_o)]
-  }
+  # stopifnot_(
+  #   checkmate::test_numeric(
+  #     x = probs, lower = 0, upper = 1, any.missing = FALSE, min.len = 1L
+  #   ),
+  #   "Argument {.arg probs} must be a {.cls numeric} vector with values
+  #     between 0 and 1."
+  # )
+  # p_i <- length(gamma_pi_raw)
+  # p_s <- length(gamma_A_raw)
+  # p_o <- length(gamma_B_raw)
+  # sds <- try(
+  #   diag(solve(-object$estimation_results$hessian)), 
+  #   silent = TRUE
+  # )
+  # if (inherits(sds, "try-error")) {
+  #   warning_(
+  #     paste0(
+  #       "Standard errors could not be computed due to singular Hessian. ",
+  #       "Confidence intervals will not be provided."
+  #     )
+  #   )
+  #   sds <- rep(NA, p_i + p_s + p_o)
+  # } else {
+  #   if (any(sds < 0)) {
+  #     warning_(
+  #       paste0(
+  #         "Standard errors could not be computed due to negative variances. ",
+  #         "Confidence intervals will not be provided."
+  #       )
+  #     )
+  #     sds <- rep(NA, p_i + p_s + p_o)
+  #   } else {
+  #     sds <- sqrt(sds)
+  #   }
+  # }
+  # for(i in seq_along(probs)) {
+  #   q <- qnorm(probs[i])
+  #   gamma_pi[paste0("q", 100 * probs[i])] <- gamma_pi_raw + q * sds[seq_len(p_i)]
+  #   gamma_A[paste0("q", 100 * probs[i])] <- gamma_A_raw + q * sds[p_i + seq_len(p_s)]
+  #   gamma_B[paste0("q", 100 * probs[i])] <- gamma_B_raw + q * sds[p_i + p_s + seq_len(p_o)]
+  # }
 
   list(
     gamma_pinitial = gamma_pi, 
@@ -101,38 +101,40 @@ coef.mnhmm <- function(object, probs = c(0.025, 0.5, 0.975), ...) {
   gamma_pi_raw <- unlist(object$coefficients$gamma_pi_raw)
   K_i <- length(object$coef_names_initial)
   gamma_pi <- data.frame(
-    state = object$state_names[-1],
+    state = unlist(lapply(object$state_names, function(x) x[-1])),
     parameter = rep(object$coef_names_initial, each = (S - 1)),
-    estimate = gamma_pi_raw,
+    estimate = unlist(gamma_pi_raw),
     cluster = rep(object$cluster_names, each = (S - 1) * K_i)
   )
   gamma_A_raw <- unlist(object$coefficients$gamma_A_raw)
   K_s <- length(object$coef_names_transition)
   gamma_A <- data.frame(
-    state_from = object$state_names,
-    state_to = rep(object$state_names[-1], each = S),
+    state_from = unlist(object$state_names),
+    state_to = rep(
+      unlist(lapply(object$state_names, function(x) x[-1])), 
+      each = S
+    ),
     parameter = rep(object$coef_names_transition, each = S * (S - 1)),
-    estimate = gamma_A_raw,
+    estimate = unlist(gamma_A_raw),
     cluster = rep(object$cluster_names, each = (S - 1) * S * K_s)
   )
-  gamma_B_raw <- unlist(object$coefficients$gamma_B_raw)
   K_o <- length(object$coef_names_emission)
   if (object$n_channels == 1) {
     gamma_B <- data.frame(
-      state = object$state_names,
+      state = unlist(object$state_names),
       observations = rep(object$symbol_names[-1], each = S),
       parameter = rep(object$coef_names_emission, each = S * (M - 1)),
-      estimate = gamma_B_raw,
+      estimate = unlist(gamma_B_raw),
       cluster =  rep(object$cluster_names, each = S * (S - 1) * K_o)
     )
   } else {
     gamma_B <- data.frame(
-      state = object$state_names,
+      state = unlist(object$state_names),
       observations = rep(unlist(lapply(object$symbol_names, "[", -1)), each = S),
       parameter = unlist(lapply(seq_len(object$n_channels), function(i) {
         rep(object$coef_names_emission, each = S * (M[i] - 1))
       })),
-      estimate = gamma_B_raw,
+      estimate = unlist(gamma_B_raw),
       cluster =  unlist(lapply(seq_len(object$n_channels), function(i) {
         rep(object$cluster_names, each = S * (M[i] - 1) * K_o)
       }))
@@ -145,36 +147,36 @@ coef.mnhmm <- function(object, probs = c(0.025, 0.5, 0.975), ...) {
     estimate = gamma_omega_raw
   )
 
-  stopifnot_(
-    checkmate::test_numeric(
-      x = probs, lower = 0, upper = 1, any.missing = FALSE, min.len = 1L
-    ),
-    "Argument {.arg probs} must be a {.cls numeric} vector with values
-      between 0 and 1."
-  )
-  p_i <- length(gamma_pi_raw)
-  p_s <- length(gamma_A_raw)
-  p_o <- length(gamma_B_raw)
-  p_d <- length(gamma_omega_raw)
-  sds <- try(
-    sqrt(diag(solve(-object$estimation_results$hessian))), 
-    silent = TRUE
-  )
-  if (inherits(sds, "try-error")) {
-    warning_(
-      "Standard errors could not be computed due to singular Hessian. 
-      Confidence intervals will not be provided."
-    )
-    sds <- rep(NA, p_i + p_s + p_o + p_d)
-  }
-
-  for(i in seq_along(probs)) {
-    q <- qnorm(probs[i])
-    gamma_pi[paste0("q", 100 * probs[i])] <- gamma_pi_raw + q * sds[seq_len(p_i)]
-    gamma_A[paste0("q", 100 * probs[i])] <- gamma_A_raw + q * sds[p_i + seq_len(p_s)]
-    gamma_B[paste0("q", 100 * probs[i])] <- gamma_B_raw + q * sds[p_i + p_s + seq_len(p_o)]
-    gamma_omega[paste0("q", 100 * probs[i])] <- gamma_omega_raw + q * sds[p_i + p_s + p_o + seq_len(p_d)]
-  }
+  # stopifnot_(
+  #   checkmate::test_numeric(
+  #     x = probs, lower = 0, upper = 1, any.missing = FALSE, min.len = 1L
+  #   ),
+  #   "Argument {.arg probs} must be a {.cls numeric} vector with values
+  #     between 0 and 1."
+  # )
+  # p_i <- length(gamma_pi_raw)
+  # p_s <- length(gamma_A_raw)
+  # p_o <- length(gamma_B_raw)
+  # p_d <- length(gamma_omega_raw)
+  # sds <- try(
+  #   sqrt(diag(solve(-object$estimation_results$hessian))), 
+  #   silent = TRUE
+  # )
+  # if (inherits(sds, "try-error")) {
+  #   warning_(
+  #     "Standard errors could not be computed due to singular Hessian. 
+  #     Confidence intervals will not be provided."
+  #   )
+  #   sds <- rep(NA, p_i + p_s + p_o + p_d)
+  # }
+  # 
+  # for(i in seq_along(probs)) {
+  #   q <- qnorm(probs[i])
+  #   gamma_pi[paste0("q", 100 * probs[i])] <- gamma_pi_raw + q * sds[seq_len(p_i)]
+  #   gamma_A[paste0("q", 100 * probs[i])] <- gamma_A_raw + q * sds[p_i + seq_len(p_s)]
+  #   gamma_B[paste0("q", 100 * probs[i])] <- gamma_B_raw + q * sds[p_i + p_s + seq_len(p_o)]
+  #   gamma_omega[paste0("q", 100 * probs[i])] <- gamma_omega_raw + q * sds[p_i + p_s + p_o + seq_len(p_d)]
+  # }
 
   list(
     gamma_initial = gamma_pi,