gamma_raw to eta, sum_to_zero parameterization of gammas, gradients t…

…o own functions
helske · Oct 3, 2024 · b5ee815 · b5ee815
1 parent 094ef02
commit b5ee815
Show file tree

Hide file tree

Showing 29 changed files with 1,475 additions and 981 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -9,40 +9,44 @@ EMx <- function(transition_, emission_, init_, obs, nSymbols, coef_, X, numberOf
     .Call(`_seqHMM_EMx`, transition_, emission_, init_, obs, nSymbols, coef_, X, numberOfStates, itermax, tol, trace, threads)
 }
 
-backward_nhmm_singlechannel <- function(gamma_A_raw, X_s, gamma_B_raw, X_o, obs) {
-    .Call(`_seqHMM_backward_nhmm_singlechannel`, gamma_A_raw, X_s, gamma_B_raw, X_o, obs)
+backward_nhmm_singlechannel <- function(eta_A, X_s, eta_B, X_o, obs) {
+    .Call(`_seqHMM_backward_nhmm_singlechannel`, eta_A, X_s, eta_B, X_o, obs)
 }
 
-backward_nhmm_multichannel <- function(gamma_A_raw, X_s, gamma_B_raw, X_o, obs, M) {
-    .Call(`_seqHMM_backward_nhmm_multichannel`, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, M)
+backward_nhmm_multichannel <- function(eta_A, X_s, eta_B, X_o, obs, M) {
+    .Call(`_seqHMM_backward_nhmm_multichannel`, eta_A, X_s, eta_B, X_o, obs, M)
 }
 
-backward_mnhmm_singlechannel <- function(gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs) {
-    .Call(`_seqHMM_backward_mnhmm_singlechannel`, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs)
+backward_mnhmm_singlechannel <- function(eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs) {
+    .Call(`_seqHMM_backward_mnhmm_singlechannel`, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs)
 }
 
-backward_mnhmm_multichannel <- function(gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs, M) {
-    .Call(`_seqHMM_backward_mnhmm_multichannel`, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs, M)
+backward_mnhmm_multichannel <- function(eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs, M) {
+    .Call(`_seqHMM_backward_mnhmm_multichannel`, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs, M)
+}
+
+cost_matrix <- function(gamma_pi_est, gamma_pi_ref, gamma_A_est, gamma_A_ref, gamma_B_est, gamma_B_ref) {
+    .Call(`_seqHMM_cost_matrix`, gamma_pi_est, gamma_pi_ref, gamma_A_est, gamma_A_ref, gamma_B_est, gamma_B_ref)
 }
 
 fast_quantiles <- function(X, probs) {
     .Call(`_seqHMM_fast_quantiles`, X, probs)
 }
 
-forward_nhmm_singlechannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs) {
-    .Call(`_seqHMM_forward_nhmm_singlechannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs)
+forward_nhmm_singlechannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs) {
+    .Call(`_seqHMM_forward_nhmm_singlechannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs)
 }
 
-forward_nhmm_multichannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, M) {
-    .Call(`_seqHMM_forward_nhmm_multichannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, M)
+forward_nhmm_multichannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs, M) {
+    .Call(`_seqHMM_forward_nhmm_multichannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs, M)
 }
 
-forward_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_forward_mnhmm_singlechannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs)
+forward_mnhmm_singlechannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs) {
+    .Call(`_seqHMM_forward_mnhmm_singlechannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs)
 }
 
-forward_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_forward_mnhmm_multichannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs, M)
+forward_mnhmm_multichannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs, M) {
+    .Call(`_seqHMM_forward_mnhmm_multichannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs, M)
 }
 
 forwardbackward <- function(transition, emission, init, obs, forwardonly, threads) {
@@ -53,36 +57,52 @@ forwardbackwardx <- function(transition, emission, init, obs, coef, X, numberOfS
     .Call(`_seqHMM_forwardbackwardx`, transition, emission, init, obs, coef, X, numberOfStates, forwardonly, threads)
 }
 
-get_omega <- function(gamma_raw, X, logspace) {
-    .Call(`_seqHMM_get_omega`, gamma_raw, X, logspace)
+get_omega <- function(gamma_raw, X) {
+    .Call(`_seqHMM_get_omega`, gamma_raw, X)
+}
+
+get_log_omega <- function(gamma_raw, X) {
+    .Call(`_seqHMM_get_log_omega`, gamma_raw, X)
+}
+
+get_omega_all <- function(gamma_raw, X) {
+    .Call(`_seqHMM_get_omega_all`, gamma_raw, X)
 }
 
-get_omega_all <- function(gamma_raw, X, logspace) {
-    .Call(`_seqHMM_get_omega_all`, gamma_raw, X, logspace)
+get_pi <- function(gamma_raw, X) {
+    .Call(`_seqHMM_get_pi`, gamma_raw, X)
 }
 
-get_pi <- function(gamma_raw, X, logspace) {
-    .Call(`_seqHMM_get_pi`, gamma_raw, X, logspace)
+get_log_pi <- function(gamma_raw, X) {
+    .Call(`_seqHMM_get_log_pi`, gamma_raw, X)
 }
 
-get_pi_all <- function(gamma_raw, X, logspace) {
-    .Call(`_seqHMM_get_pi_all`, gamma_raw, X, logspace)
+get_pi_all <- function(gamma_raw, X) {
+    .Call(`_seqHMM_get_pi_all`, gamma_raw, X)
 }
 
-get_A <- function(gamma_raw, X, logspace, tv) {
-    .Call(`_seqHMM_get_A`, gamma_raw, X, logspace, tv)
+get_A <- function(gamma_raw, X, tv) {
+    .Call(`_seqHMM_get_A`, gamma_raw, X, tv)
 }
 
-get_A_all <- function(gamma_raw, X, logspace, tv) {
-    .Call(`_seqHMM_get_A_all`, gamma_raw, X, logspace, tv)
+get_log_A <- function(gamma_raw, X, tv) {
+    .Call(`_seqHMM_get_log_A`, gamma_raw, X, tv)
 }
 
-get_B <- function(gamma_raw, X, logspace, add_missing, tv) {
-    .Call(`_seqHMM_get_B`, gamma_raw, X, logspace, add_missing, tv)
+get_A_all <- function(gamma_raw, X, tv) {
+    .Call(`_seqHMM_get_A_all`, gamma_raw, X, tv)
 }
 
-get_B_all <- function(gamma_raw, X, logspace, add_missing, tv) {
-    .Call(`_seqHMM_get_B_all`, gamma_raw, X, logspace, add_missing, tv)
+get_B <- function(gamma_raw, X, add_missing, tv) {
+    .Call(`_seqHMM_get_B`, gamma_raw, X, add_missing, tv)
+}
+
+get_log_B <- function(gamma_raw, X, add_missing, tv) {
+    .Call(`_seqHMM_get_log_B`, gamma_raw, X, add_missing, tv)
+}
+
+get_B_all <- function(gamma_raw, X, add_missing, tv) {
+    .Call(`_seqHMM_get_B_all`, gamma_raw, X, add_missing, tv)
 }
 
 logLikHMM <- function(transition, emission, init, obs, threads) {
@@ -125,20 +145,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, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti) {
-    .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, Ti)
+log_objective_nhmm_singlechannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti) {
+    .Call(`_seqHMM_log_objective_nhmm_singlechannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti)
 }
 
-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, Ti) {
-    .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, Ti)
+log_objective_nhmm_multichannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs, M, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti) {
+    .Call(`_seqHMM_log_objective_nhmm_multichannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs, M, iv_pi, iv_A, iv_B, tv_A, tv_B, Ti)
 }
 
-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, Ti) {
-    .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, Ti)
+log_objective_mnhmm_singlechannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs, iv_pi, iv_A, iv_B, tv_A, tv_B, iv_omega, Ti) {
+    .Call(`_seqHMM_log_objective_mnhmm_singlechannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs, iv_pi, iv_A, iv_B, tv_A, tv_B, iv_omega, Ti)
 }
 
-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, Ti) {
-    .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, Ti)
+log_objective_mnhmm_multichannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs, M, iv_pi, iv_A, iv_B, tv_A, tv_B, iv_omega, Ti) {
+    .Call(`_seqHMM_log_objective_mnhmm_multichannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs, M, iv_pi, iv_A, iv_B, tv_A, tv_B, iv_omega, Ti)
 }
 
 log_objectivex <- function(transition, emission, init, obs, ANZ, BNZ, INZ, nSymbols, coef, X, numberOfStates, threads) {
@@ -153,8 +173,12 @@ objectivex <- function(transition, emission, init, obs, ANZ, BNZ, INZ, nSymbols,
     .Call(`_seqHMM_objectivex`, transition, emission, init, obs, ANZ, BNZ, INZ, nSymbols, coef, X, numberOfStates, threads)
 }
 
-softmax <- function(x, logspace) {
-    .Call(`_seqHMM_softmax`, x, logspace)
+softmax <- function(x) {
+    .Call(`_seqHMM_softmax`, x)
+}
+
+eta_to_gamma <- function(x) {
+    .Call(`_seqHMM_eta_to_gamma`, x)
 }
 
 varcoef <- function(coef, X) {
@@ -165,20 +189,20 @@ viterbi <- function(transition, emission, init, obs) {
     .Call(`_seqHMM_viterbi`, transition, emission, init, obs)
 }
 
-viterbi_nhmm_singlechannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs) {
-    .Call(`_seqHMM_viterbi_nhmm_singlechannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs)
+viterbi_nhmm_singlechannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs) {
+    .Call(`_seqHMM_viterbi_nhmm_singlechannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs)
 }
 
-viterbi_nhmm_multichannel <- function(gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, M) {
-    .Call(`_seqHMM_viterbi_nhmm_multichannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, obs, M)
+viterbi_nhmm_multichannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs, M) {
+    .Call(`_seqHMM_viterbi_nhmm_multichannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, obs, M)
 }
 
-viterbi_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_viterbi_mnhmm_singlechannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs)
+viterbi_mnhmm_singlechannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs) {
+    .Call(`_seqHMM_viterbi_mnhmm_singlechannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs)
 }
 
-viterbi_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_viterbi_mnhmm_multichannel`, gamma_pi_raw, X_i, gamma_A_raw, X_s, gamma_B_raw, X_o, gamma_omega_raw, X_d, obs, M)
+viterbi_mnhmm_multichannel <- function(eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs, M) {
+    .Call(`_seqHMM_viterbi_mnhmm_multichannel`, eta_pi, X_i, eta_A, X_s, eta_B, X_o, eta_omega, X_d, obs, M)
 }
 
 viterbix <- function(transition, emission, init, obs, coef, X, numberOfStates) {

diff --git a/R/bootstrap.R b/R/bootstrap.R
@@ -17,19 +17,29 @@ bootstrap_model <- function(model) {
 }
 #' @export
 bootstrap_coefs.nhmm <- function(model, B = 1000, 
-                                 method = c("nonparametric", "parametric")) {
+                                 method = c("nonparametric", "parametric"),
+                                 penalty) {
   method <- match.arg(method)
   stopifnot_(
     checkmate::test_int(x = B, lower = 0L), 
     "Argument {.arg B} must be a single positive integer."
   )
   init <- model$coefficients
-  penalty <- model$estimation_results$penalty
+  if (missing(penalty)) {
+    penalty <- model$estimation_results$penalty
+  }
+  mle_coefficients <- model$coefficients
   if (method == "nonparametric") {
     coefs <- matrix(NA, length(unlist(init)), B)
     for (i in seq_len(B)) {
       mod <- bootstrap_model(model)
       fit <- fit_nhmm(mod, init, 0, 0, 1, penalty, FALSE)
+      m <- cost_matrix(fit$coefficients$eta_pi, mle_coefficients$eta_pi,
+                       fit$coefficients$eta_A, mle_coefficients$eta_A,
+                       fit$coefficients$eta_B, mle_coefficients$eta_B,
+                       fit$X_initial, fit$X_transition, fit$X_emission)
+      perm <- RcppHungarian:HungarianSolver(m)$pairs[, 2]
+      fit$coefficients$gamma_pi[perm]
       coefs[, i] <- unlist(fit$coefficients)
     }
   } else {
@@ -63,7 +73,9 @@ bootstrap_coefs.mnhmm <- function(model, B = 1000,
     "Argument {.arg B} must be a single positive integer."
   )
   init <- model$coefficients
-  penalty <- model$estimation_results$penalty
+  if (missing(penalty)) {
+    penalty <- model$estimation_results$penalty
+  }
   if (method == "nonparametric") {
     coefs <- matrix(NA, length(unlist(init)), B)
     for (i in seq_len(B)) {