add bounds, remove pseudocounts

helske · Dec 3, 2024 · 2240195 · 2240195
1 parent 72cbbe2
commit 2240195
Show file tree

Hide file tree

Showing 24 changed files with 167 additions and 274 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -129,20 +129,20 @@ logSumExp <- function(x) {
     .Call(`_seqHMM_logSumExp`, x)
 }
 
-EM_LBFGS_mnhmm_singlechannel <- function(eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, pseudocount, bound) {
-    .Call(`_seqHMM_EM_LBFGS_mnhmm_singlechannel`, eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, pseudocount, bound)
+EM_LBFGS_mnhmm_singlechannel <- function(eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) {
+    .Call(`_seqHMM_EM_LBFGS_mnhmm_singlechannel`, eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound)
 }
 
-EM_LBFGS_mnhmm_multichannel <- function(eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, pseudocount, bound) {
-    .Call(`_seqHMM_EM_LBFGS_mnhmm_multichannel`, eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, pseudocount, bound)
+EM_LBFGS_mnhmm_multichannel <- function(eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) {
+    .Call(`_seqHMM_EM_LBFGS_mnhmm_multichannel`, eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound)
 }
 
-EM_LBFGS_nhmm_singlechannel <- function(eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, pseudocount, bound) {
-    .Call(`_seqHMM_EM_LBFGS_nhmm_singlechannel`, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, pseudocount, bound)
+EM_LBFGS_nhmm_singlechannel <- function(eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) {
+    .Call(`_seqHMM_EM_LBFGS_nhmm_singlechannel`, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound)
 }
 
-EM_LBFGS_nhmm_multichannel <- function(eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, pseudocount, bound) {
-    .Call(`_seqHMM_EM_LBFGS_nhmm_multichannel`, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, pseudocount, bound)
+EM_LBFGS_nhmm_multichannel <- function(eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) {
+    .Call(`_seqHMM_EM_LBFGS_nhmm_multichannel`, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound)
 }
 
 backward_nhmm_singlechannel <- function(eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B) {

diff --git a/R/bootstrap.R b/R/bootstrap.R
@@ -99,7 +99,6 @@ bootstrap_coefs.nhmm <- function(model, nsim = 1000,
   init <- model$etas
   gammas_mle <- model$gammas
   lambda <- model$estimation_results$lambda
-  pseudocount <- model$estimation_results$pseudocount
   bound <- model$estimation_results$bound
   p <- progressr::progressor(along = seq_len(nsim))
   original_options <- options(future.globals.maxSize = Inf)
@@ -109,8 +108,7 @@ bootstrap_coefs.nhmm <- function(model, nsim = 1000,
       seq_len(nsim), function(i) {
         mod <- bootstrap_model(model)
         fit <- fit_nhmm(mod, init, init_sd = 0, restarts = 0, lambda = lambda, 
-                        method = method, pseudocount = pseudocount, 
-                        bound = bound, ...)
+                        method = method, bound = bound, ...)
         if (fit$estimation_results$return_code >= 0) {
           fit$gammas <- permute_states(fit$gammas, gammas_mle)
         } else {
@@ -137,8 +135,7 @@ bootstrap_coefs.nhmm <- function(model, nsim = 1000,
           N, T_, M, S, formula_pi, formula_A, formula_B,
           d, time, id, init, 0)$model
         fit <- fit_nhmm(mod, init, init_sd = 0, restarts = 0, lambda = lambda, 
-                        method = method, pseudocount = pseudocount, 
-                        bound = bound, ...)
+                        method = method, bound = bound, ...)
         if (fit$estimation_results$return_code >= 0) {
           fit$gammas <- permute_states(fit$gammas, gammas_mle)
         } else {
@@ -185,7 +182,6 @@ bootstrap_coefs.mnhmm <- function(model, nsim = 1000,
   gammas_mle <- model$gammas
   pcp_mle <- posterior_cluster_probabilities(model)
   lambda <- model$estimation_results$lambda
-  pseudocount <- model$estimation_results$pseudocount
   bound <- model$estimation_results$bound
   D <- model$n_clusters
   p <- progressr::progressor(along = seq_len(nsim))
@@ -196,8 +192,7 @@ bootstrap_coefs.mnhmm <- function(model, nsim = 1000,
       seq_len(nsim), function(i) {
         mod <- bootstrap_model(model)
         fit <- fit_mnhmm(mod, init, init_sd = 0, restarts = 0, lambda = lambda, 
-                         method = method, pseudocount = pseudocount, 
-                         bound = bound, ...)
+                         method = method, bound = bound, ...)
         if (fit$estimation_results$return_code >= 0) {
           fit <- permute_clusters(fit, pcp_mle)
           for (j in seq_len(D)) {
@@ -234,8 +229,7 @@ bootstrap_coefs.mnhmm <- function(model, nsim = 1000,
           N, T_, M, S, D, formula_pi, formula_A, formula_B, formula_omega,
           d, time, id, init, 0)$model
         fit <- fit_mnhmm(mod, init, init_sd = 0, restarts = 0, lambda = lambda, 
-                         method = method, pseudocount = pseudocount, 
-                         bound = bound, ...)
+                         method = method, bound = bound, ...)
         if (fit$estimation_results$return_code >= 0) {
           fit <- permute_clusters(fit, pcp_mle)
           for (j in seq_len(D)) {

diff --git a/R/dnm_mnhmm.R b/R/dnm_mnhmm.R
@@ -132,8 +132,8 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
     out <- future.apply::future_lapply(seq_len(restarts), function(i) {
       init <- unlist(create_initial_values(inits, model, init_sd))
       fit <- nloptr(
-        x0 = init, eval_f = objectivef, lb = -bound, ub = bound,
-        opts = control_restart
+        x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
+        ub = rep(bound, length(init)), opts = control_restart
       )
       p()
       fit
@@ -162,8 +162,8 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
   }
 
   out <- nloptr(
-    x0 = init, eval_f = objectivef, lb = -bound, ub = bound,
-    opts = control
+    x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
+    ub = rep(bound, length(init)), opts = control
   )
   if (out$status < 0) {
     warning_(
@@ -205,7 +205,6 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
     return_codes_of_restarts = if(restarts > 0L) return_codes else NULL,
     all_solutions = all_solutions,
     lambda = lambda,
-    pseudocount = 0,
     bound = bound,
     method = "DNM", 
     algorithm = control$algorithm

diff --git a/R/dnm_nhmm.R b/R/dnm_nhmm.R
@@ -98,8 +98,8 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
     out <- future.apply::future_lapply(seq_len(restarts), function(i) {
       init <- unlist(create_initial_values(inits, model, init_sd))
       fit <- nloptr(
-        x0 = init, eval_f = objectivef, lb = -bound, ub = bound,
-        opts = control_restart
+        x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
+        ub = rep(bound, length(init)), opts = control_restart
       )
       p()
       fit
@@ -128,8 +128,8 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
   }
 
   out <- nloptr(
-    x0 = init, eval_f = objectivef, lb = -bound, ub = bound, 
-    opts = control
+    x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
+    ub = rep(bound, length(init)), opts = control
   )
   if (out$status < 0) {
     warning_(
@@ -165,7 +165,6 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
     return_codes_of_restarts = if(restarts > 0L) return_codes else NULL,
     all_solutions = all_solutions,
     lambda = lambda,
-    pseudocount = 0,
     bound = bound,
     method = "DNM", 
     algorithm = control$algorithm

diff --git a/R/em_dnm_mnhmm.R b/R/em_dnm_mnhmm.R
@@ -1,4 +1,4 @@
-em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount, 
+em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, 
                          bound, control, control_restart, control_mstep, 
                          save_all_solutions) {
   M <- model$n_symbols
@@ -143,7 +143,7 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
           control_restart$print_level, control_mstep$maxeval,
           control_mstep$ftol_abs, control_mstep$ftol_rel,
           control_mstep$xtol_abs, control_mstep$xtol_rel, 
-          control_mstep$print_level, lambda, pseudocount, bound)
+          control_mstep$print_level, lambda, bound)
       } else {
         eta_B <- unlist(init$B, recursive = FALSE)
         fit <- EM_LBFGS_mnhmm_multichannel(
@@ -156,7 +156,7 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
           control_restart$print_level, control_mstep$maxeval,
           control_mstep$ftol_abs, control_mstep$ftol_rel,
           control_mstep$xtol_abs, control_mstep$xtol_rel, 
-          control_mstep$print_level, lambda, pseudocount, bound)
+          control_mstep$print_level, lambda, bound)
       }
       em_return_code <- fit$return_code
       if (em_return_code >= 0) {
@@ -171,8 +171,8 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
           )
         )
         fit <- nloptr(
-          x0 = init, eval_f = objectivef, 
-          lb = -bound, ub = bound, opts = control_restart
+          x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
+          ub = rep(bound, length(init)), opts = control_restart
         )
         p()
         fit
@@ -229,7 +229,7 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
         control$print_level, control_mstep$maxeval,
         control_mstep$ftol_abs, control_mstep$ftol_rel,
         control_mstep$xtol_abs, control_mstep$xtol_rel, 
-        control_mstep$print_level, lambda, pseudocount, bound)
+        control_mstep$print_level, lambda, bound)
     } else {
       eta_B <- unlist(init$B, recursive = FALSE)
       out <- EM_LBFGS_mnhmm_multichannel(
@@ -242,7 +242,7 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
         control$print_level, control_mstep$maxeval,
         control_mstep$ftol_abs, control_mstep$ftol_rel,
         control_mstep$xtol_abs, control_mstep$xtol_rel, 
-        control_mstep$print_level, lambda, pseudocount, bound)
+        control_mstep$print_level, lambda, bound)
     }
     em_return_code <- out$return_code
     if (em_return_code >= 0) {
@@ -263,8 +263,9 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
     }
   }
   out <- nloptr(
-    x0 = unlist(init), eval_f = objectivef, lb = -bound, ub = bound,
-    opts = control
+    x0 = unlist(init), eval_f = objectivef, 
+    lb = -rep(bound, length(unlist(init))), 
+    ub = rep(bound, length(unlist(init))), opts = control
   )
   if (out$status < 0) {
     warning_(
@@ -306,7 +307,6 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
     return_codes_of_restarts = if(restarts > 0L) return_codes else NULL,
     all_solutions = all_solutions,
     lambda = lambda,
-    pseudocount = pseudocount,
     bound = bound,
     method = "EM-DNM",
     algorithm = control$algorithm,

diff --git a/R/em_dnm_nhmm.R b/R/em_dnm_nhmm.R
@@ -1,4 +1,4 @@
-em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount, 
+em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda,
                         bound, control, control_restart, control_mstep, 
                         save_all_solutions) {
   M <- model$n_symbols
@@ -108,7 +108,7 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
           control_restart$print_level, control_mstep$maxeval,
           control_mstep$ftol_abs, control_mstep$ftol_rel,
           control_mstep$xtol_abs, control_mstep$xtol_rel, 
-          control_mstep$print_level, lambda, pseudocount, bound)
+          control_mstep$print_level, lambda, bound)
       } else {
         fit <- EM_LBFGS_nhmm_multichannel(
           init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B, obs,
@@ -119,7 +119,7 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
           control_restart$print_level, control_mstep$maxeval,
           control_mstep$ftol_abs, control_mstep$ftol_rel,
           control_mstep$xtol_abs, control_mstep$xtol_rel, 
-          control_mstep$print_level, lambda, pseudocount, bound)
+          control_mstep$print_level, lambda, bound)
       }
       em_return_code <- fit$return_code
       if (em_return_code >= 0) {
@@ -133,8 +133,8 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
           )
         )
         fit <- nloptr(
-          x0 = init, eval_f = objectivef, lb = -bound, ub = bound, 
-          opts = control_restart
+          x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
+          ub = rep(bound, length(init)), opts = control_restart
         )
         p()
         fit
@@ -187,7 +187,7 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
         control$print_level, control_mstep$maxeval,
         control_mstep$ftol_abs, control_mstep$ftol_rel,
         control_mstep$xtol_abs, control_mstep$xtol_rel, 
-        control_mstep$print_level, lambda, pseudocount, bound)
+        control_mstep$print_level, lambda, bound)
     } else {
       out <- EM_LBFGS_nhmm_multichannel(
         init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B, obs,
@@ -198,7 +198,7 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
         control$print_level, control_mstep$maxeval,
         control_mstep$ftol_abs, control_mstep$ftol_rel,
         control_mstep$xtol_abs, control_mstep$xtol_rel, 
-        control_mstep$print_level, lambda, pseudocount, bound)
+        control_mstep$print_level, lambda, bound)
     }
     em_return_code <- out$return_code
     if (em_return_code >= 0) {
@@ -218,8 +218,9 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
     }
   }
   out <- nloptr(
-    x0 = unlist(init), eval_f = objectivef, lb = -bound, ub = bound,
-    opts = control
+    x0 = unlist(init), eval_f = objectivef, 
+    lb = -rep(bound, length(unlist(init))), 
+    ub = rep(bound, length(unlist(init))), opts = control
   )
   if (out$status < 0) {
     warning_(
@@ -252,7 +253,6 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
     return_codes_of_restarts = if(restarts > 0L) return_codes else NULL,
     all_solutions = all_solutions,
     lambda = lambda,
-    pseudocount = pseudocount,
     bound = bound,
     method = "EM-DNM",
     algorithm = control$algorithm,