fix initial values based on probs

helske · Oct 5, 2024 · 741b8bc · 741b8bc
1 parent 440f3d8
commit 741b8bc
Show file tree

Hide file tree

Showing 4 changed files with 31 additions and 30 deletions.
diff --git a/R/create_initial_values.R b/R/create_initial_values.R
@@ -185,14 +185,14 @@ create_eta_omega_inits <- function(x, D, K, init_sd = 0) {
 #' @noRd
 create_inits_vector <- function(x, n, K, sd = 0, D = 1) {
   cbind(
-    inv_softmax(x)[-1], 
+    p_to_eta(x), # intercepts
     matrix(rnorm((n - 1) * (K - 1), sd = sd), n - 1, K - 1)
   )
 }
 create_inits_matrix <- function(x, n, m, K, sd = 0) {
   z <- array(0, c(m - 1, K, n))
   for (i in seq_len(n)) {
-    z[, , i] <- t(create_inits_vector(x[i, ], m, K, sd))
+    z[, , i] <- create_inits_vector(x[i, ], m, K, sd)
   }
   z
 }

diff --git a/R/fit_mnhmm.R b/R/fit_mnhmm.R
@@ -91,9 +91,9 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
     Qd <- t(create_Q(D))
     if (need_grad) {
       objectivef <- function(pars) {
-        if (any(!is.finite(exp(pars)))) {
-          return(list(objective = Inf, gradient = rep(-Inf, length(pars))))
-        } 
+        # if (any(!is.finite(exp(pars)))) {
+        #   return(list(objective = Inf, gradient = rep(-Inf, length(pars))))
+        # } 
         eta_pi <- create_eta_pi_mnhmm(pars[seq_len(n_i)], S, K_i, D)
         eta_A <- create_eta_A_mnhmm(pars[n_i + seq_len(n_s)], S, K_s, D)
         eta_B <- create_eta_B_mnhmm(
@@ -112,9 +112,9 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
       }
     } else {
       objectivef <- function(pars) {
-        if (any(!is.finite(exp(pars)))) {
-          return(Inf)
-        } 
+        # if (any(!is.finite(exp(pars)))) {
+        #   return(Inf)
+        # } 
         eta_pi <- create_eta_pi_mnhmm(pars[seq_len(n_i)], S, K_i, D)
         eta_A <- create_eta_A_mnhmm(pars[n_i + seq_len(n_s)], S, K_s, D)
         eta_B <- create_eta_B_mnhmm(
@@ -137,9 +137,9 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
     Qd <- t(create_Q(D))
     if (need_grad) {
       objectivef <- function(pars) {
-        if (any(!is.finite(exp(pars)))) {
-          return(list(objective = Inf, gradient = rep(-Inf, length(pars))))
-        }
+        # if (any(!is.finite(exp(pars)))) {
+        #   return(list(objective = Inf, gradient = rep(-Inf, length(pars))))
+        # }
         eta_pi <- create_eta_pi_mnhmm(pars[seq_len(n_i)], S, K_i, D)
         eta_A <- create_eta_A_mnhmm(
           pars[n_i + seq_len(n_s)], 
@@ -164,9 +164,9 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
       }
     } else {
       objectivef <- function(pars) {
-        if (any(!is.finite(exp(pars)))) {
-          return(Inf)
-        } 
+        # if (any(!is.finite(exp(pars)))) {
+        #   return(Inf)
+        # } 
         eta_pi <- create_eta_pi_mnhmm(pars[seq_len(n_i)], S, K_i, D)
         eta_A <- create_eta_A_mnhmm(
           pars[n_i + seq_len(n_s)], S, K_s, D

diff --git a/R/fit_nhmm.R b/R/fit_nhmm.R
@@ -77,9 +77,9 @@ fit_nhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
     Qm <- t(create_Q(M))
     if (need_grad) {
       objectivef <- function(pars) {
-        if (any(!is.finite(exp(pars)))) {
-          return(list(objective = Inf, gradient = rep(-Inf, length(pars))))
-        }
+        # if (any(!is.finite(exp(pars)))) {
+        #   return(list(objective = Inf, gradient = rep(-Inf, length(pars))))
+        # }
         eta_pi <- create_eta_pi_nhmm(pars[seq_len(n_i)], S, K_i)
         eta_A <- create_eta_A_nhmm(pars[n_i + seq_len(n_s)], S, K_s)
         eta_B <- create_eta_B_nhmm(
@@ -94,9 +94,9 @@ fit_nhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
       }
     } else {
       objectivef <- function(pars) {
-        if (any(!is.finite(exp(pars)))) {
-          return(Inf)
-        } 
+        # if (any(!is.finite(exp(pars)))) {
+        #   return(Inf)
+        # } 
         eta_pi <- create_eta_pi_nhmm(pars[seq_len(n_i)], S, K_i)
         eta_A <- create_eta_A_nhmm(pars[n_i + seq_len(n_s)], S, K_s)
         eta_B <- create_eta_B_nhmm(
@@ -114,9 +114,9 @@ fit_nhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
     Qm <- lapply(M, function(m) t(create_Q(m)))
     if (need_grad) {
       objectivef <- function(pars) {
-        if (any(!is.finite(exp(pars)))) {
-          return(list(objective = Inf, gradient = rep(-Inf, length(pars))))
-        }
+        # if (any(!is.finite(exp(pars)))) {
+        #   return(list(objective = Inf, gradient = rep(-Inf, length(pars))))
+        # }
         eta_pi <- create_eta_pi_nhmm(pars[seq_len(n_i)], S, K_i)
         eta_A <- create_eta_A_nhmm(pars[n_i + seq_len(n_s)], S, K_s)
         eta_B <- create_eta_multichannel_B_nhmm(
@@ -131,9 +131,9 @@ fit_nhmm <- function(model, inits, init_sd, restarts, threads, penalty, ...) {
       }
     } else {
       objectivef <- function(pars) {
-        if (any(!is.finite(exp(pars)))) {
-          return(Inf)
-        }
+        # if (any(!is.finite(exp(pars)))) {
+        #   return(Inf)
+        # }
         eta_pi <- create_eta_pi_nhmm(pars[seq_len(n_i)], S, K_i)
         eta_A <- create_eta_A_nhmm(pars[n_i + seq_len(n_s)], S, K_s)
         eta_B <- create_eta_multichannel_B_nhmm(

diff --git a/R/utilities.R b/R/utilities.R
@@ -17,15 +17,16 @@ is_list_of_lists <- function(x) {
     }
   }
 }
-#' Regularized Inverse Softmax Function
+#' (Regularized) Inverse of softmax(Q*eta)
 #' 
 #' @noRd
-inv_softmax <- function(x) {
+p_to_eta <-function(x) {
+  Q <- create_Q(length(x))
   x <- pmin(pmax(x, 0.001), 0.999)
   x <- x / sum(x)
-  log(x) - log(x[1])
+  log_x <- log(x)
+  t(Q) %*% (log_x - mean(log_x))
 }
-
 #' Stop Function Execution Unless Condition Is True
 #' 
 #' Function copied from the `dynamite` package.