diff --git a/R/bootstrap.R b/R/bootstrap.R
index f6bd68f..403dfb9 100644
--- a/R/bootstrap.R
+++ b/R/bootstrap.R
@@ -76,7 +76,7 @@ permute_clusters <- function(model, pcp_mle) {
 #' sequences with replacement, whereas the latter simulates new datasets based 
 #' on the model.
 #' @param method Estimation method used in bootstrapping. Defaults to 
-#' `"EM-LBFGS"`.
+#' `"EM-DNM"`.
 #' @param ... Additional arguments to [estimate_nhmm()] or [estimate_mnhmm()].
 #' @return The original model with additional element `model$boot`, which 
 #' contains A n * B matrix where n is the number of coefficients. The order of 
@@ -90,7 +90,7 @@ bootstrap_coefs <- function(model, ...) {
 #' @export
 bootstrap_coefs.nhmm <- function(model, B = 1000, 
                                  type = c("nonparametric", "parametric"),
-                                 method = "EM-LBFGS", ...) {
+                                 method = "EM-DNM", ...) {
   type <- match.arg(type)
   stopifnot_(
     checkmate::test_int(x = B, lower = 0L), 
@@ -147,7 +147,7 @@ bootstrap_coefs.nhmm <- function(model, B = 1000,
 #' @export
 bootstrap_coefs.mnhmm <- function(model, B = 1000, 
                                   type = c("nonparametric", "parametric"),
-                                  method = "EM-LBFGS", ...) {
+                                  method = "EM-DNM", ...) {
   type <- match.arg(type)
   stopifnot_(
     checkmate::test_int(x = B, lower = 0L), 
diff --git a/R/dnm_mnhmm.R b/R/dnm_mnhmm.R
index ba78443..ae7710b 100644
--- a/R/dnm_mnhmm.R
+++ b/R/dnm_mnhmm.R
@@ -1,4 +1,4 @@
-dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, control, 
+dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, control, 
                        control_restart, save_all_solutions) {
   M <- model$n_symbols
   S <- model$n_states
diff --git a/R/em_lbfgs_mnhmm.R b/R/em_dnm_mnhmm.R
similarity index 91%
rename from R/em_lbfgs_mnhmm.R
rename to R/em_dnm_mnhmm.R
index c9b9c38..dd01eb5 100644
--- a/R/em_lbfgs_mnhmm.R
+++ b/R/em_dnm_mnhmm.R
@@ -1,6 +1,6 @@
-em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount, 
-                          control, control_restart, control_mstep, 
-                          save_all_solutions) {
+em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount, 
+                           control, control_restart, control_mstep, 
+                           save_all_solutions) {
   M <- model$n_symbols
   S <- model$n_states
   T_ <- model$length_of_sequences
@@ -191,7 +191,22 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
       )
     }
     optimum <- successful[which.max(logliks[successful])]
-    init <- out[[optimum]]$solution
+    pars <- out[[optimum]]$solution
+    init <- list()
+    init$pi <- create_eta_pi_mnhmm(pars[seq_len(n_i)], S, K_pi, D)
+    init$A <- create_eta_A_mnhmm(pars[n_i + seq_len(n_s)], S, K_A, D)
+    if (C == 1L) {
+      init$B <- create_eta_B_mnhmm(
+        pars[n_i + n_s + seq_len(n_o)], S, M, K_B, D
+      )
+    } else {
+      init$B <- create_eta_multichannel_B_mnhmm(
+        pars[n_i + n_s + seq_len(n_o)], S, M, K_B, D
+      )
+    }
+    init$omega <- create_eta_omega_mnhmm(
+      pars[n_i + n_s + n_o + seq_len(n_d)], D, K_omega
+    )
     if (save_all_solutions) {
       all_solutions <- out
     }
@@ -228,7 +243,7 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
     init <- unlist(
       create_initial_values(
         stats::setNames(
-          fit[c("eta_omega", "eta_pi", "eta_A", "eta_B")], 
+          out[c("eta_omega", "eta_pi", "eta_A", "eta_B")], 
           c("omega", "pi", "A", "B")
         ), 
         model, 
@@ -238,7 +253,7 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
   } else {
     warning_(
       paste("EM-step terminated due to error:", error_msg(out$return_code),
-            "Running L-BFGS using initial values for EM.")
+            "Running DNM using initial values for EM.")
     )
     init <- unlist(create_initial_values(inits, model, init_sd))
   }
@@ -286,7 +301,8 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
     time = end_time - start_time,
     lambda = lambda,
     pseudocount = 0,
-    method = "EM-LBFGS"
+    method = "EM-DNM",
+    algorithm = control$algorithm
   )
   model
-}
\ No newline at end of file
+}
diff --git a/R/em_lbfgs_nhmm.R b/R/em_dnm_nhmm.R
similarity index 91%
rename from R/em_lbfgs_nhmm.R
rename to R/em_dnm_nhmm.R
index 53a4f27..755735e 100644
--- a/R/em_lbfgs_nhmm.R
+++ b/R/em_dnm_nhmm.R
@@ -1,4 +1,4 @@
-em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount, 
+em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount, 
                           control, control_restart, control_mstep, 
                           save_all_solutions) {
   M <- model$n_symbols
@@ -102,7 +102,7 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
         fit <- EM_LBFGS_nhmm_singlechannel(
           init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B, obs,
           Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B,
-          n_obs, control$maxeval_em_lbfgs,
+          n_obs, control$maxeval_em_dnm,
           control_restart$ftol_abs, control_restart$ftol_rel,
           control_restart$xtol_abs, control_restart$xtol_rel, 
           control_restart$print_level, control_mstep$maxeval,
@@ -113,7 +113,7 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
         fit <- EM_LBFGS_nhmm_multichannel(
           init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B, obs,
           Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B,
-          n_obs, control$maxeval_em_lbfgs,
+          n_obs, control$maxeval_em_dnm,
           control_restart$ftol_abs, control_restart$ftol_rel,
           control_restart$xtol_abs, control_restart$xtol_rel, 
           control_restart$print_level, control_mstep$maxeval,
@@ -153,7 +153,19 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
       )
     }
     optimum <- successful[which.max(logliks[successful])]
-    init <- out[[optimum]]$solution
+    pars <- out[[optimum]]$solution
+    init <- list()
+    init$pi <- create_eta_pi_nhmm(pars[seq_len(n_i)], S, K_pi)
+    init$A <- create_eta_A_nhmm(pars[n_i + seq_len(n_s)], S, K_A)
+    if (C == 1L) {
+      init$B <- create_eta_B_nhmm(
+        pars[n_i + n_s + seq_len(n_o)], S, M, K_B
+      )
+    } else {
+      init$B <- create_eta_multichannel_B_nhmm(
+        pars[n_i + n_s + seq_len(n_o)], S, M, K_B
+      )
+    }
     if (save_all_solutions) {
       all_solutions <- out
     }
@@ -164,7 +176,7 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
     out <- EM_LBFGS_nhmm_singlechannel(
       init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B, obs,
       Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B,
-      n_obs, control$maxeval_em_lbfgs,
+      n_obs, control$maxeval_em_dnm,
       control$ftol_abs, control$ftol_rel,
       control$xtol_abs, control$xtol_rel, 
       control$print_level, control_mstep$maxeval,
@@ -175,7 +187,7 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
     out <- EM_LBFGS_nhmm_multichannel(
       init$pi, model$X_pi, init$A, model$X_A, init$B, model$X_B, obs,
       Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B,
-      n_obs, control$maxeval_em_lbfgs,
+      n_obs, control$maxeval_em_dnm,
       control$ftol_abs, control$ftol_rel,
       control$xtol_abs, control$xtol_rel, 
       control$print_level, control_mstep$maxeval,
@@ -187,7 +199,7 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
     init <- unlist(
       create_initial_values(
         stats::setNames(
-          fit[c("eta_pi", "eta_A", "eta_B")], c("pi", "A", "B")
+          out[c("eta_pi", "eta_A", "eta_B")], c("pi", "A", "B")
         ), 
         model, 
         init_sd = 0
@@ -196,7 +208,7 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
   } else {
     warning_(
       paste("EM-step terminated due to error:", error_msg(out$return_code),
-            "Running L-BFGS using initial values for EM.")
+            "Running DNM using initial values for EM.")
     )
     init <- unlist(create_initial_values(inits, model, init_sd))
   }
@@ -238,7 +250,8 @@ em_lbfgs_nhmm <- function(model, inits, init_sd, restarts, lambda, pseudocount,
     time = end_time - start_time,
     lambda = lambda,
     pseudocount = 0,
-    method = "EM-LBFGS"
+    method = "EM-DNM",
+    algorithm = control$algorithm
   )
   model
-}
\ No newline at end of file
+}
diff --git a/R/estimate_mnhmm.R b/R/estimate_mnhmm.R
index 9c31e46..708de6c 100644
--- a/R/estimate_mnhmm.R
+++ b/R/estimate_mnhmm.R
@@ -44,11 +44,10 @@ 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 = 0, method = "EM-LBFGS", pseudocount = 1e-4, 
+    restarts = 0L, lambda = 0, method = "EM-DNM", pseudocount = 1e-4, 
     store_data = TRUE, ...) {
   
   call <- match.call()
-  method <- match.arg(method, c("EM-LBFGS", "DNM", "EM"))
   model <- build_mnhmm(
     observations, n_states, n_clusters, initial_formula, 
     transition_formula, emission_formula, cluster_formula, data, time, id, 
@@ -58,10 +57,6 @@ estimate_mnhmm <- function(
     checkmate::test_flag(x = store_data), 
     "Argument {.arg store_data} must be a single {.cls logical} value."
   )
-  stopifnot_(
-    checkmate::check_number(lambda, lower = 0), 
-    "Argument {.arg lambda} must be a single non-negative {.cls numeric} value."
-  )
   if (store_data) {
     model$data <- data
   }
diff --git a/R/estimate_nhmm.R b/R/estimate_nhmm.R
index ff4742c..473e7ad 100644
--- a/R/estimate_nhmm.R
+++ b/R/estimate_nhmm.R
@@ -4,7 +4,7 @@
 #' `nhmm` where initial, transition and emission probabilities 
 #' (potentially) depend on covariates.
 #' 
-#' By default, the model parameters are estimated using EM-LBFGS algorithm 
+#' By default, the model parameters are estimated using EM-DNM algorithm 
 #' which first runs some iterations (100 by default) of EM algorithm, and then 
 #' switches to L-BFGS. Other options include any numerical optimization 
 #' algorithm of [nloptr::nloptr()], or plain EM algorithm where the 
@@ -28,14 +28,14 @@
 #' `...`.  These, as well as argument `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-LBFGS` 
-#' can be set using argument `maxeval_em_lbfgs` (default is 100), and 
+#' 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). 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 sepately for the M-step of EM algorithm via list 
+#' options can be defined separately for the M-step of EM algorithm via list 
 #' `control_em`. 
 #' 
 #' @param observations Either the name of the response variable in `data`, or 
@@ -78,8 +78,8 @@
 #' @param method Optimization method used. Option `"EM"` uses EM
 #' algorithm with L-BFGS in the M-step. Option `"DNM"` uses 
 #' direct maximization of the log-likelihood, by default using L-BFGS. Option 
-#' `"EM-LBFGS"` (the default) runs first a maximum of 100 iterations of EM and 
-#' then switches to L-BFGS.
+#' `"EM-DNM"` (the default) runs first a maximum of 100 iterations of EM and 
+#' then switches to L-BFGS (but other algorithms of NLopt can be used).
 #' @param pseudocount A positive scalar to be added for the expected counts of 
 #' E-step. Only used in EM algorithm. Default is 1e-4. Larger values can be used 
 #' to avoid zero probabilities in initial, transition, and emission 
@@ -111,11 +111,10 @@ 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 = 0, method = "EM-LBFGS", pseudocount = 1e-4, store_data = TRUE, 
+    lambda = 0, method = "EM-DNM", pseudocount = 1e-4, store_data = TRUE, 
     ...) {
   
   call <- match.call()
-  method <- match.arg(method, c("EM-LBFGS", "DNM", "EM"))
   model <- build_nhmm(
     observations, n_states, initial_formula, 
     transition_formula, emission_formula, data, time, id, state_names, 
@@ -125,10 +124,6 @@ estimate_nhmm <- function(
     checkmate::test_flag(x = store_data), 
     "Argument {.arg store_data} must be a single {.cls logical} value."
   )
-  stopifnot_(
-    checkmate::check_number(lambda, lower = 0), 
-    "Argument {.arg lambda} must be a single non-negative {.cls numeric} value."
-  )
   if (store_data) {
     model$data <- data
   }
diff --git a/R/fit_mnhmm.R b/R/fit_mnhmm.R
index 699b7a5..6c9dc8a 100644
--- a/R/fit_mnhmm.R
+++ b/R/fit_mnhmm.R
@@ -9,6 +9,11 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method,
     checkmate::test_int(x = restarts, lower = 0L), 
     "Argument {.arg restarts} must be a single integer."
   )
+  method <- match.arg(method, c("EM-DNM", "DNM", "EM"))
+  stopifnot_(
+    checkmate::check_number(lambda, lower = 0), 
+    "Argument {.arg lambda} must be a single non-negative {.cls numeric} value."
+  )
   control <- utils::modifyList(
     list(
       ftol_abs = 1e-8,
@@ -18,7 +23,7 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method,
       maxeval = 1e4,
       print_level = 0,
       algorithm = "NLOPT_LD_LBFGS",
-      maxeval_em_lbfgs = 100
+      maxeval_em_dnm = 100
     ),
     list(...)
   )
@@ -57,10 +62,10 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method,
     return(model)
   }
   all_solutions <- NULL
-  if (method == "EM-LBFGS") {
-    out <- lbfgs_em_mnhmm(
-      model, inits, init_sd, restarts, lambda, control, control_restart, 
-      save_all_solutions 
+  if (method == "EM-DNM") {
+    out <- em_dnm_mnhmm(
+      model, inits, init_sd, restarts, lambda, pseudocount, control, 
+      control_restart, control_mstep, save_all_solutions
     )
   }
   if (method == "DNM") {
diff --git a/R/fit_nhmm.R b/R/fit_nhmm.R
index a62fea7..375bb05 100644
--- a/R/fit_nhmm.R
+++ b/R/fit_nhmm.R
@@ -9,6 +9,11 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method, pseudocoun
     checkmate::test_int(x = restarts, lower = 0L), 
     "Argument {.arg restarts} must be a single integer."
   )
+  method <- match.arg(method, c("EM-DNM", "DNM", "EM"))
+  stopifnot_(
+    checkmate::check_number(lambda, lower = 0), 
+    "Argument {.arg lambda} must be a single non-negative {.cls numeric} value."
+  )
   control <- utils::modifyList(
     list(
       ftol_abs = 1e-8,
@@ -18,7 +23,7 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method, pseudocoun
       maxeval = 1e4,
       print_level = 0,
       algorithm = "NLOPT_LD_LBFGS",
-      maxeval_em_lbfgs = 100
+      maxeval_em_dnm = 100
     ),
     list(...)
   )
@@ -53,8 +58,8 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method, pseudocoun
     }
     return(model)
   }
-  if (method == "EM-LBFGS") {
-    out <- lbfgs_em_nhmm(
+  if (method == "EM-DNM") {
+    out <- em_dnm_nhmm(
       model, inits, init_sd, restarts, lambda, pseudocount, control, 
       control_restart, control_mstep, save_all_solutions
     )
diff --git a/man/bootstrap.Rd b/man/bootstrap.Rd
index 674d36d..c051052 100644
--- a/man/bootstrap.Rd
+++ b/man/bootstrap.Rd
@@ -12,7 +12,7 @@ bootstrap_coefs(model, ...)
   model,
   B = 1000,
   type = c("nonparametric", "parametric"),
-  method = "EM-LBFGS",
+  method = "EM-DNM",
   ...
 )
 
@@ -20,7 +20,7 @@ bootstrap_coefs(model, ...)
   model,
   B = 1000,
   type = c("nonparametric", "parametric"),
-  method = "EM-LBFGS",
+  method = "EM-DNM",
   ...
 )
 }
@@ -37,7 +37,7 @@ sequences with replacement, whereas the latter simulates new datasets based
 on the model.}
 
 \item{method}{Estimation method used in bootstrapping. Defaults to
-\code{"EM-LBFGS"}.}
+\code{"EM-DNM"}.}
 }
 \value{
 The original model with additional element \code{model$boot}, which
diff --git a/man/estimate_mnhmm.Rd b/man/estimate_mnhmm.Rd
index 7c99a07..353e63d 100644
--- a/man/estimate_mnhmm.Rd
+++ b/man/estimate_mnhmm.Rd
@@ -22,7 +22,7 @@ estimate_mnhmm(
   init_sd = 2,
   restarts = 0L,
   lambda = 0,
-  method = "EM-LBFGS",
+  method = "EM-DNM",
   pseudocount = 1e-04,
   store_data = TRUE,
   ...
@@ -93,8 +93,8 @@ the penalization term is scaled with number of non-missing observations.}
 \item{method}{Optimization method used. Option \code{"EM"} uses EM
 algorithm with L-BFGS in the M-step. Option \code{"DNM"} uses
 direct maximization of the log-likelihood, by default using L-BFGS. Option
-\code{"EM-LBFGS"} (the default) runs first a maximum of 100 iterations of EM and
-then switches to L-BFGS.}
+\code{"EM-DNM"} (the default) runs first a maximum of 100 iterations of EM and
+then switches to L-BFGS (but other algorithms of NLopt can be used).}
 
 \item{pseudocount}{A positive scalar to be added for the expected counts of
 E-step. Only used in EM algorithm. Default is 1e-4. Larger values can be used
diff --git a/man/estimate_nhmm.Rd b/man/estimate_nhmm.Rd
index f0ec1f2..fd8ec10 100644
--- a/man/estimate_nhmm.Rd
+++ b/man/estimate_nhmm.Rd
@@ -19,7 +19,7 @@ estimate_nhmm(
   init_sd = 2,
   restarts = 0L,
   lambda = 0,
-  method = "EM-LBFGS",
+  method = "EM-DNM",
   pseudocount = 1e-04,
   store_data = TRUE,
   ...
@@ -80,8 +80,8 @@ the penalization term is scaled with number of non-missing observations.}
 \item{method}{Optimization method used. Option \code{"EM"} uses EM
 algorithm with L-BFGS in the M-step. Option \code{"DNM"} uses
 direct maximization of the log-likelihood, by default using L-BFGS. Option
-\code{"EM-LBFGS"} (the default) runs first a maximum of 100 iterations of EM and
-then switches to L-BFGS.}
+\code{"EM-DNM"} (the default) runs first a maximum of 100 iterations of EM and
+then switches to L-BFGS (but other algorithms of NLopt can be used).}
 
 \item{pseudocount}{A positive scalar to be added for the expected counts of
 E-step. Only used in EM algorithm. Default is 1e-4. Larger values can be used
@@ -105,7 +105,7 @@ Function \code{estimate_nhmm} estimates a hidden Markov model object of class
 (potentially) depend on covariates.
 }
 \details{
-By default, the model parameters are estimated using EM-LBFGS algorithm
+By default, the model parameters are estimated using EM-DNM algorithm
 which first runs some iterations (100 by default) of EM algorithm, and then
 switches to L-BFGS. Other options include any numerical optimization
 algorithm of \code{\link[nloptr:nloptr]{nloptr::nloptr()}}, or plain EM algorithm where the
@@ -129,14 +129,14 @@ by passing arguments \code{ftol_abs}, \code{ftol_rel}, \code{xtol_abs}, and \cod
 \code{...}.  These, as well as argument \code{maxeval} (maximum number of iterations,
 1e4 by default), and \code{print_level} (default is \code{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 \code{EM-LBFGS}
-can be set using argument \code{maxeval_em_lbfgs} (default is 100), and
+main optimization method. The number of initial EM iterations in \code{EM-DNM}
+can be set using argument \code{maxeval_em_dnm} (default is 100), and
 algorithm for direct numerical optimization can be defined using argument
 \code{algorithm} (see \code{\link[nloptr:nloptr]{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 \code{control_restart} which takes
 a list such as \code{list(ftol_rel = 0.01, print_level = 1)}. Additionally, same
-options can be defined sepately for the M-step of EM algorithm via list
+options can be defined separately for the M-step of EM algorithm via list
 \code{control_em}.
 }
 \examples{
diff --git a/src/mnhmm_EM.cpp b/src/mnhmm_EM.cpp
index d4d1aff..71fcdca 100644
--- a/src/mnhmm_EM.cpp
+++ b/src/mnhmm_EM.cpp
@@ -86,7 +86,7 @@ void mnhmm_base::mstep_omega(const double xtol_abs, const double ftol_abs,
   double minf;
   mstep_iter = 0;
   int status = nlopt_optimize(opt_omega, x_omega.memptr(), &minf);
-  if (print_level > 0 && status > 0) {
+  if (print_level > 2 && status > 0) {
     Rcpp::Rcout<<"M-step of cluster probabilities ended with status "<<status<<
       " after "<<mstep_iter<<" iterations."<<std::endl;
   }
@@ -182,7 +182,7 @@ void mnhmm_base::mstep_pi(const double xtol_abs, const double ftol_abs,
     x_pi = arma::vectorise(eta_pi(d));
     mstep_iter = 0;
     status = nlopt_optimize(opt_pi, x_pi.memptr(), &minf);
-    if (print_level > 0 && status > 0) {
+    if (print_level > 2 && status > 0) {
       Rcpp::Rcout<<"M-step of initial probabilities ended with status "<<status<<
         " after "<<mstep_iter<<" iterations."<<std::endl;
     }
@@ -302,7 +302,7 @@ void mnhmm_base::mstep_A(const double ftol_abs, const double ftol_rel,
       x_A = arma::vectorise(eta_A(d).slice(s));
       mstep_iter = 0;
       status = nlopt_optimize(opt_A, x_A.memptr(), &minf);
-      if (print_level > 0 && status > 0) {
+      if (print_level > 2 && status > 0) {
         Rcpp::Rcout<<"M-step of transition probabilities of state "<<s + 1<<
           " ended with status "<<status<<" after "<<mstep_iter<<
             " iterations."<<std::endl;
@@ -429,7 +429,7 @@ void mnhmm_sc::mstep_B(const double ftol_abs, const double ftol_rel,
       x_B = arma::vectorise(eta_B(d).slice(s));
       mstep_iter = 0;
       status = nlopt_optimize(opt_B, x_B.memptr(), &minf);
-      if (print_level > 0 && status > 0) {
+      if (print_level > 2 && status > 0) {
         Rcpp::Rcout<<"M-step of emission probabilities of state "<<s + 1<<
           " ended with status "<<status<<" after "<<mstep_iter<<
             " iterations."<<std::endl;
@@ -560,7 +560,7 @@ void mnhmm_mc::mstep_B(const double ftol_abs, const double ftol_rel,
         x_B = arma::vectorise(eta_B(c, d).slice(s));
         mstep_iter = 0;
         status = nlopt_optimize(opt_B, x_B.memptr(), &minf);
-        if (print_level > 0 && status > 0) {
+        if (print_level > 2 && status > 0) {
           Rcpp::Rcout<<"M-step of emission probabilities of state "<<s + 1<<
             " and channel "<<c<<" ended with status "<<status<<" after "<<mstep_iter<<
               " iterations."<<std::endl;
diff --git a/src/nhmm_EM.cpp b/src/nhmm_EM.cpp
index 933c6fd..7131a7c 100644
--- a/src/nhmm_EM.cpp
+++ b/src/nhmm_EM.cpp
@@ -86,7 +86,7 @@ void nhmm_base::mstep_pi(const double xtol_abs, const double ftol_abs,
   double minf;
   mstep_iter = 0;
   int status = nlopt_optimize(opt_pi, x_pi.memptr(), &minf);
-  if (print_level > 0 && status > 0) {
+  if (print_level > 2 && status > 0) {
     Rcpp::Rcout<<"M-step of initial probabilities ended with status "<<status<<
       " after "<<mstep_iter<<" iterations."<<std::endl;
   }
@@ -202,7 +202,7 @@ void nhmm_base::mstep_A(const double ftol_abs, const double ftol_rel,
     x_A = arma::vectorise(eta_A.slice(s));
     mstep_iter = 0;
     status = nlopt_optimize(opt_A, x_A.memptr(), &minf);
-    if (print_level > 0 && status > 0) {
+    if (print_level > 2 && status > 0) {
       Rcpp::Rcout<<"M-step of transition probabilities of state "<<s + 1<<
         " ended with status "<<status<<" after "<<mstep_iter<<
           " iterations."<<std::endl;
@@ -323,7 +323,7 @@ void nhmm_sc::mstep_B(const double ftol_abs, const double ftol_rel,
     x_B = arma::vectorise(eta_B.slice(s));
     mstep_iter = 0;
     status = nlopt_optimize(opt_B, x_B.memptr(), &minf);
-    if (print_level > 0 && status > 0) {
+    if (print_level > 2 && status > 0) {
       Rcpp::Rcout<<"M-step of emission probabilities of state "<<s + 1<<
         " ended with status "<<status<<" after "<<mstep_iter<<
           " iterations."<<std::endl;
@@ -450,7 +450,7 @@ void nhmm_mc::mstep_B(const double ftol_abs, const double ftol_rel,
       x_B = arma::vectorise(eta_B(c).slice(s));
       mstep_iter = 0;
       status = nlopt_optimize(opt_B, x_B.memptr(), &minf);
-      if (print_level > 0 && status > 0) {
+      if (print_level > 2 && status > 0) {
         Rcpp::Rcout<<"M-step of emission probabilities of state "<<s + 1<<
           " and channel "<<c<<" ended with status "<<status<<" after "<<mstep_iter<<
             " iterations."<<std::endl;