check for gradient norm

R/RcppExports.R

@@ -137,12 +137,12 @@ EM_LBFGS_mnhmm_multichannel <- function(eta_omega, X_omega, eta_pi, X_pi, eta_A,
     .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, 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_singlechannel <- function(eta_pi, x, 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, 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, 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)
+EM_LBFGS_nhmm_multichannel <- function(eta_pi, x, 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, 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) {

R/dnm_mnhmm.R

@@ -135,6 +135,10 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
         x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
         ub = rep(bound, length(init)), opts = control_restart
       )
+      if (fit$status == -1 && need_grad) {
+        grad_norm <- sqrt(sum(objectivef(fit$solution)$gradient^2))
+        if (grad_norm < 1e-6) fit$status <- 6
+      }
       p()
       fit
     },
@@ -145,8 +149,8 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
     successful <- which(return_codes > 0)
     if (length(successful) == 0) {
       warning_(
-        c("All optimizations terminated due to error.",
-          "Error of first restart: ", error_msg(return_codes[1]),
+        c("All restarts terminated due to error.",
+          "Error of first restart: ", return_msg(return_codes[1]),
           "Trying once more. ")
       )
       init <- unlist(create_initial_values(inits, model, init_sd))
@@ -165,9 +169,13 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
     x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
     ub = rep(bound, length(init)), opts = control
   )
+  if (ou$status == -1 && need_grad) {
+    grad_norm <- sqrt(sum(objectivef(out$solution)$gradient^2))
+    if (grad_norm < 1e-6) out$status <- 6
+  }
   if (out$status < 0) {
     warning_(
-      paste("Optimization terminated due to error:", error_msg(out$status))
+      paste("Optimization terminated due to error:", return_msg(out$status))
     )
     loglik <- NaN
   } else {
@@ -210,4 +218,4 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
     algorithm = control$algorithm
   )
   model
-}
+}

R/dnm_nhmm.R

@@ -101,6 +101,10 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
         x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
         ub = rep(bound, length(init)), opts = control_restart
       )
+      if (fit$status == -1 && need_grad) {
+        grad_norm <- sqrt(sum(objectivef(fit$solution)$gradient^2))
+        if (grad_norm < 1e-6) fit$status <- 6
+      }
       p()
       fit
     },
@@ -111,8 +115,8 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
     successful <- which(return_codes > 0)
     if (length(successful) == 0) {
       warning_(
-        c("All optimizations terminated due to error.",
-          "Error of first restart: ", error_msg(return_codes[1]),
+        c("All restarts terminated due to error.",
+          "Error of first restart: ", return_msg(return_codes[1]),
           "Trying once more. ")
       )
       init <- unlist(create_initial_values(inits, model, init_sd))
@@ -131,9 +135,13 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
     x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
     ub = rep(bound, length(init)), opts = control
   )
+  if (out$status == -1 && need_grad) {
+    grad_norm <- sqrt(sum(objectivef(out$solution)$gradient^2))
+    if (grad_norm < 1e-6) out$status <- 6
+  }
   if (out$status < 0) {
     warning_(
-      paste("Optimization terminated due to error:", error_msg(out$status))
+      paste("Optimization terminated due to error:", return_msg(out$status))
     )
     loglik <- NaN
   } else {
@@ -170,4 +178,4 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control,
     algorithm = control$algorithm
   )
   model
-}
+}

R/em_dnm_mnhmm.R

@@ -174,6 +174,10 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda,
           x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
           ub = rep(bound, length(init)), opts = control_restart
         )
+        if (fit$status == -1 && need_grad) {
+          grad_norm <- sqrt(sum(objectivef(fit$solution)$gradient^2))
+          if (grad_norm < 1e-6) fit$status <- 6
+        }
         p()
         fit
       } else {
@@ -187,8 +191,8 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda,
     successful <- which(return_codes > 0)
     if (length(successful) == 0) {
       warning_(
-        c("All optimizations terminated due to error.",
-          "Error of first restart: ", error_msg(return_codes[1]),
+        c("All restarts terminated due to error.",
+          "Error of first restart: ", return_msg(return_codes[1]),
           "Running DNM using initial values for EM.")
       )
       init <- create_initial_values(inits, model, init_sd)
@@ -256,7 +260,7 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda,
       )
     } else {
       warning_(
-        paste("EM-step terminated due to error:", error_msg(em_return_code),
+        paste("EM-step terminated due to error:", return_msg(em_return_code),
               "Running DNM using initial values for EM.")
       )
       init <- create_initial_values(inits, model, init_sd)
@@ -267,9 +271,13 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda,
     lb = -rep(bound, length(unlist(init))), 
     ub = rep(bound, length(unlist(init))), opts = control
   )
+  if (out$status == -1 && need_grad) {
+    grad_norm <- sqrt(sum(objectivef(out$solution)$gradient^2))
+    if (grad_norm < 1e-6) out$status <- 6
+  }
   if (out$status < 0) {
     warning_(
-      paste("Optimization terminated due to error:", error_msg(out$status))
+      paste("Optimization terminated due to error:", return_msg(out$status))
     )
     loglik <- NaN
   } else {

R/em_dnm_nhmm.R

@@ -136,6 +136,10 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda,
           x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), 
           ub = rep(bound, length(init)), opts = control_restart
         )
+        if (fit$status == -1 && need_grad) {
+          grad_norm <- sqrt(sum(objectivef(fit$solution)$gradient^2))
+          if (grad_norm < 1e-6) fit$status <- 6
+        }
         p()
         fit
       } else {
@@ -149,8 +153,8 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda,
     successful <- which(return_codes > 0)
     if (length(successful) == 0) {
       warning_(
-        c("All optimizations terminated due to error.",
-          "Error of first restart: ", error_msg(return_codes[1]),
+        c("All restarts terminated due to error.",
+          "Error of first restart: ", return_msg(return_codes[1]),
           "Running DNM using initial values for EM.")
       )
       em_return_code <- return_codes[1] + 1000
@@ -211,7 +215,7 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda,
       )
     } else {
       warning_(
-        paste("EM-step terminated due to error:", error_msg(em_return_code),
+        paste("EM-step terminated due to error:", return_msg(em_return_code),
               "Running DNM using initial values for EM.")
       )
       init <- create_initial_values(inits, model, init_sd)
@@ -222,9 +226,13 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda,
     lb = -rep(bound, length(unlist(init))), 
     ub = rep(bound, length(unlist(init))), opts = control
   )
+  if (out$status == -1 && need_grad) {
+    grad_norm <- sqrt(sum(objectivef(out$solution)$gradient^2))
+    if (grad_norm < 1e-6) out$status <- 6
+  }
   if (out$status < 0) {
     warning_(
-      paste("Optimization terminated due to error:", error_msg(out$status))
+      paste("Optimization terminated due to error:", return_msg(out$status))
     )
   }
 

R/em_mnhmm.R

@@ -72,8 +72,8 @@ em_mnhmm <- function(model, inits, init_sd, restarts, lambda,
     return_codes <- unlist(lapply(out, "[[", "return_code"))
     if (all(return_codes < 0)) {
       warning_(
-        c("All optimizations terminated due to error.",
-          "Error of first restart: ", error_msg(return_codes[1]))
+        c("All restarts terminated due to error.",
+          "Error of first restart: ", return_msg(return_codes[1]))
       )
     }
     logliks <- unlist(lapply(out, "[[", "logLik"))
@@ -116,7 +116,7 @@ em_mnhmm <- function(model, inits, init_sd, restarts, lambda,
   }
   if (out$return_code < 0) {
     warning_(
-      paste("Optimization terminated due to error:", error_msg(out$return_code))
+      paste("Optimization terminated due to error:", return_msg(out$return_code))
     )
   }
   model$etas$pi <- c(out$eta_pi)
@@ -152,4 +152,4 @@ em_mnhmm <- function(model, inits, init_sd, restarts, lambda,
     method = "EM"
   )
   model
-}
+}

R/em_nhmm.R

@@ -64,8 +64,8 @@ em_nhmm <- function(model, inits, init_sd, restarts, lambda,
     return_codes <- unlist(lapply(out, "[[", "return_code"))
     if (all(return_codes < 0)) {
       warning_(
-        c("All optimizations terminated due to error.",
-          "Error of first restart: ", error_msg(return_codes[1]))
+        c("All restarts terminated due to error.",
+          "Error of first restart: ", return_msg(return_codes[1]))
       )
     }
     logliks <- unlist(lapply(out, "[[", "logLik"))
@@ -104,7 +104,7 @@ em_nhmm <- function(model, inits, init_sd, restarts, lambda,
   }
   if (out$return_code < 0) {
     warning_(
-      paste("Optimization terminated due to error:", error_msg(out$return_code))
+      paste("Optimization terminated due to error:", return_msg(out$return_code))
     )
   }
   model$etas$pi[] <- out$eta_pi
@@ -136,4 +136,4 @@ em_nhmm <- function(model, inits, init_sd, restarts, lambda,
     method = "EM"
   )
   model
-}
+}

R/estimate_mnhmm.R

@@ -44,8 +44,8 @@ estimate_mnhmm <- function(
     transition_formula = ~1, emission_formula = ~1, cluster_formula = ~1,
     data = NULL, time = NULL, id = NULL, state_names = NULL, 
     channel_names = NULL, cluster_names = NULL, inits = "random", init_sd = 2, 
-    restarts = 0L, lambda = 1e-4, method = "EM-DNM", 
-    bound = 50, store_data = TRUE, ...) {
+    restarts = 0L, lambda = 0, method = "EM-DNM", bound = 50, 
+    control_restart = list(), control_mstep = list(), store_data = TRUE, ...) {
 
   call <- match.call()
   model <- build_mnhmm(
@@ -60,8 +60,10 @@ estimate_mnhmm <- function(
   if (store_data) {
     model$data <- data
   }
+  control <- list(...)
   start_time <- proc.time()
-  out <- fit_mnhmm(model, inits, init_sd, restarts, lambda, method, bound, ...)
+  out <- fit_mnhmm(model, inits, init_sd, restarts, lambda, method, bound, 
+                   control, control_restart, control_mstep)
   end_time <- proc.time()
   out$estimation_results$time <- end_time - start_time
   attr(out, "call") <- call

R/estimate_nhmm.R

@@ -22,22 +22,26 @@
 #' standardardized before optimization.
 #'  
 #' By default, the convergence is claimed when the relative 
-#' change of the objective function is less than `1e-10`, or the 
+#' change of the objective function is less than `1e-8`, or the 
 #' relative change of the parameters is less than `1e-6`. These can be changed
 #' by passing arguments `ftol_rel` and `xtol_rel` via `...`. These, as well as 
-#' arguments `ftol_abs` and `xtol_abs`  for absolute changes, `maxeval` 
-#' (maximum number of iterations, 1e4 by default), and `print_level` (default 
-#' is `0`, no console output of  optimization, larger values are more verbose), 
-#' are used by the chosen main optimization method. The number of initial EM 
-#' iterations in `EM-DNM` can be set using argument `maxeval_em_dnm` (default 
-#' is 10), and algorithm for direct numerical optimization can be defined using 
-#' argument `algorithm` (see [nloptr::nloptr()] for possible options). It is also 
-#' possible to separately control these stopping criteria for the multistart 
-#' phase by defining (some of) them via argument `control_restart` which takes 
-#' a list such as `list(ftol_rel = 0.01, print_level = 1)`. Additionally, same 
-#' options can be defined separately for the M-step of EM algorithm via list 
-#' `control_mstep`. By default, `control_restart` and `control_mstep` match the 
-#' the main options defined via `...`. 
+#' arguments `ftol_abs` and `xtol_abs` for absolute changes 
+#' (`0` by default), `maxeval` (maximum number of iterations, 1e4 by default), 
+#' and `print_level` (default is `0`, no console output of  optimization, 
+#' larger values are more verbose), are used by the chosen main optimization 
+#' method. The number of initial EM iterations in `EM-DNM` can be set using 
+#' argument `maxeval_em_dnm` (default is 100), and algorithm for direct
+#' numerical optimization can be defined using argument `algorithm` 
+#' (see [nloptr::nloptr()] for possible options).
+#' 
+#' For controlling these stopping criteria for the multistart phase, argument 
+#' `control_restart` takes a list such as `list(ftol_rel = 0.01, print_level = 1)`. 
+#' Default are as in the case of main optimization (which is always run once a
+#' fter the restarts, using best solution from restarts as initial value),
+#' except that `ftol_rel = 1e-6` and that `xtol_rel = 1e-4`.
+#' Additionally, same options can be defined separately for the M-step of EM 
+#' algorithm via list `control_mstep`. For `control_mstep`, the 
+#' default value are `ftol_rel = 1e-6`, `xtol_rel = 1e-4`, and `maxeval = 100`.
 #' 
 #' @references Steven G. Johnson, The NLopt nonlinear-optimization package, http://github.com/stevengj/nlopt
 #' @param observations Either the name of the response variable in `data`, or 
@@ -117,9 +121,8 @@ estimate_nhmm <- function(
     transition_formula = ~1, emission_formula = ~1, 
     data = NULL, time = NULL, id = NULL, state_names = NULL, 
     channel_names = NULL, inits = "random", init_sd = 2, restarts = 0L, 
-    lambda = 1e-4, method = "EM-DNM", bound = 50, 
-    store_data = TRUE, 
-    ...) {
+    lambda = 0, method = "EM-DNM", bound = 50, control_restart = list(), 
+    control_mstep = list(), store_data = TRUE, ...) {
 
   call <- match.call()
   model <- build_nhmm(
@@ -134,8 +137,10 @@ estimate_nhmm <- function(
   if (store_data) {
     model$data <- data
   }
+  control <- list(...)
   start_time <- proc.time()
-  out <- fit_nhmm(model, inits, init_sd, restarts, lambda, method, bound, ...)
+  out <- fit_nhmm(model, inits, init_sd, restarts, lambda, method, bound, 
+                  control, control_restart, control_mstep)
   end_time <- proc.time()
   out$estimation_results$time <- end_time - start_time
   attr(out, "call") <- call