Skip to content

Commit

Permalink
new sum-to-zero parameterization, modularise get_probs etc
Browse files Browse the repository at this point in the history
  • Loading branch information
@helske
helske committed Oct 4, 2024
1 parent b5ee815 commit be42ffa
Show file tree
Hide file tree
Showing 39 changed files with 950 additions and 753 deletions.
2 changes: 2 additions & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,8 @@ export("cluster_names<-")
export("state_names<-")
export(alphabet)
export(ame)
export(bootstrap_coefs.mnhmm)
export(bootstrap_coefs.nhmm)
export(build_hmm)
export(build_lcm)
export(build_mhmm)
Expand Down
96 changes: 56 additions & 40 deletions R/RcppExports.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,18 +17,38 @@ 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(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_singlechannel <- function(eta_A, X_s, eta_B, X_o, obs) {
.Call(`_seqHMM_backward_mnhmm_singlechannel`, eta_A, X_s, eta_B, X_o, obs)
}

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)
backward_mnhmm_multichannel <- function(eta_A, X_s, eta_B, X_o, obs, M) {
.Call(`_seqHMM_backward_mnhmm_multichannel`, eta_A, X_s, eta_B, X_o, 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)
}

create_Q <- function(n) {
.Call(`_seqHMM_create_Q`, n)
}

eta_to_gamma_mat <- function(eta) {
.Call(`_seqHMM_eta_to_gamma_mat`, eta)
}

eta_to_gamma_cube <- function(eta) {
.Call(`_seqHMM_eta_to_gamma_cube`, eta)
}

eta_to_gamma_mat_field <- function(eta) {
.Call(`_seqHMM_eta_to_gamma_mat_field`, eta)
}

eta_to_gamma_cube_field <- function(eta) {
.Call(`_seqHMM_eta_to_gamma_cube_field`, eta)
}

fast_quantiles <- function(X, probs) {
.Call(`_seqHMM_fast_quantiles`, X, probs)
}
Expand Down Expand Up @@ -57,52 +77,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) {
.Call(`_seqHMM_get_omega`, gamma_raw, X)
get_omega <- function(gamma, X) {
.Call(`_seqHMM_get_omega`, gamma, X)
}

get_log_omega <- function(gamma_raw, X) {
.Call(`_seqHMM_get_log_omega`, gamma_raw, X)
get_log_omega <- function(gamma, X) {
.Call(`_seqHMM_get_log_omega`, gamma, X)
}

get_omega_all <- function(gamma_raw, X) {
.Call(`_seqHMM_get_omega_all`, gamma_raw, X)
get_omega_all <- function(gamma, X) {
.Call(`_seqHMM_get_omega_all`, gamma, X)
}

get_pi <- function(gamma_raw, X) {
.Call(`_seqHMM_get_pi`, gamma_raw, X)
get_pi <- function(gamma, X) {
.Call(`_seqHMM_get_pi`, gamma, X)
}

get_log_pi <- function(gamma_raw, X) {
.Call(`_seqHMM_get_log_pi`, gamma_raw, X)
get_log_pi <- function(gamma, X) {
.Call(`_seqHMM_get_log_pi`, gamma, X)
}

get_pi_all <- function(gamma_raw, X) {
.Call(`_seqHMM_get_pi_all`, gamma_raw, X)
get_pi_all <- function(gamma, X) {
.Call(`_seqHMM_get_pi_all`, gamma, X)
}

get_A <- function(gamma_raw, X, tv) {
.Call(`_seqHMM_get_A`, gamma_raw, X, tv)
get_A <- function(gamma, X, tv) {
.Call(`_seqHMM_get_A`, gamma, X, tv)
}

get_log_A <- function(gamma_raw, X, tv) {
.Call(`_seqHMM_get_log_A`, gamma_raw, X, tv)
get_log_A <- function(gamma, X, tv) {
.Call(`_seqHMM_get_log_A`, gamma, X, tv)
}

get_A_all <- function(gamma_raw, X, tv) {
.Call(`_seqHMM_get_A_all`, gamma_raw, X, tv)
get_A_all <- function(gamma, X, tv) {
.Call(`_seqHMM_get_A_all`, gamma, X, tv)
}

get_B <- function(gamma_raw, X, add_missing, tv) {
.Call(`_seqHMM_get_B`, gamma_raw, X, add_missing, tv)
get_B <- function(gamma, X, add_missing, tv) {
.Call(`_seqHMM_get_B`, gamma, 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_log_B <- function(gamma, X, add_missing, tv) {
.Call(`_seqHMM_get_log_B`, gamma, 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)
get_B_all <- function(gamma, X, add_missing, tv) {
.Call(`_seqHMM_get_B_all`, gamma, X, add_missing, tv)
}

logLikHMM <- function(transition, emission, init, obs, threads) {
Expand Down Expand Up @@ -145,20 +165,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(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_singlechannel <- function(Qs, Qm, 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`, Qs, Qm, 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(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_nhmm_multichannel <- function(Qs, Qm, 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`, Qs, Qm, 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(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_singlechannel <- function(Qs, Qm, Qd, 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`, Qs, Qm, Qd, 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(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_objective_mnhmm_multichannel <- function(Qs, Qm, Qd, 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`, Qs, Qm, Qd, 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) {
Expand All @@ -177,10 +197,6 @@ softmax <- function(x) {
.Call(`_seqHMM_softmax`, x)
}

eta_to_gamma <- function(x) {
.Call(`_seqHMM_eta_to_gamma`, x)
}

varcoef <- function(coef, X) {
.Call(`_seqHMM_varcoef`, coef, X)
}
Expand Down
110 changes: 50 additions & 60 deletions R/coef.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,36 +10,32 @@
coef.nhmm <- function(object, probs = c(0.025, 0.975), ...) {
S <- object$n_states
M <- object$n_symbols
gamma_pi <- c(eta_to_gamma(object$coefficients$eta_pi))
gamma_pi <- data.frame(
state = object$state_names,
parameter = rep(object$coef_names_initial, each = S),
estimate = gamma_pi
estimate = c(object$gammas$pi)
)
eta_A <- c(object$coefficients$eta_A)
gamma_A <- data.frame(
state_from = object$state_names,
state_to = rep(object$state_names[-1], each = S),
parameter = rep(object$coef_names_transition, each = S * (S - 1)),
estimate = eta_A
state_to = rep(object$state_names, each = S),
parameter = rep(object$coef_names_transition, each = S^2),
estimate = c(object$gammas$A)
)
if (object$n_channels == 1) {
eta_B <- c(object$coefficients$eta_B)
gamma_B <- data.frame(
state = object$state_names,
observation = rep(object$symbol_names[-1], each = S),
parameter = rep(object$coef_names_emission, each = S * (M - 1)),
estimate = eta_B
observation = rep(object$symbol_names, each = S),
parameter = rep(object$coef_names_emission, each = S * M),
estimate = c(object$gammas$B)
)
} else {
eta_B <- unlist(object$coefficients$eta_B)
gamma_B <- data.frame(
state = object$state_names,
observation = rep(unlist(lapply(object$symbol_names, "[", -1)), each = S),
observation = rep(unlist(object$symbol_names), each = S),
parameter = unlist(lapply(seq_len(object$n_channels), function(i) {
rep(object$coef_names_emission, each = S * (M[i] - 1))
rep(object$coef_names_emission, each = S * M[i])
})),
estimate = eta_B
estimate = unlist(object$gammas$B)
)
}

Expand All @@ -50,9 +46,9 @@ coef.nhmm <- function(object, probs = c(0.025, 0.975), ...) {
# "Argument {.arg probs} must be a {.cls numeric} vector with values
# between 0 and 1."
# )
# p_i <- length(eta_pi)
# p_s <- length(eta_A)
# p_o <- length(eta_B)
# p_i <- length(gamma_pi)
# p_s <- length(gamma_A)
# p_o <- length(gamma_B)
# sds <- try(
# diag(solve(-object$estimation_results$hessian)),
# silent = TRUE
Expand Down Expand Up @@ -80,15 +76,15 @@ coef.nhmm <- function(object, probs = c(0.025, 0.975), ...) {
# }
# for(i in seq_along(probs)) {
# q <- qnorm(probs[i])
# gamma_pi[paste0("q", 100 * probs[i])] <- eta_pi + q * sds[seq_len(p_i)]
# gamma_A[paste0("q", 100 * probs[i])] <- eta_A + q * sds[p_i + seq_len(p_s)]
# gamma_B[paste0("q", 100 * probs[i])] <- eta_B + q * sds[p_i + p_s + seq_len(p_o)]
# gamma_pi[paste0("q", 100 * probs[i])] <- gamma_pi + q * sds[seq_len(p_i)]
# gamma_A[paste0("q", 100 * probs[i])] <- gamma_A + q * sds[p_i + seq_len(p_s)]
# gamma_B[paste0("q", 100 * probs[i])] <- gamma_B + q * sds[p_i + p_s + seq_len(p_o)]
# }

list(
gamma_initial = gamma_pi,
beta_transition = gamma_A,
beta_emission = gamma_B
initial = gamma_pi,
transition = gamma_A,
emission = gamma_B
)
}
#' @rdname coef
Expand All @@ -98,54 +94,48 @@ coef.mnhmm <- function(object, probs = c(0.025, 0.5, 0.975), ...) {
S <- object$n_states
M <- object$n_symbols
D <- object$n_clusters
eta_pi <- unlist(object$coefficients$eta_pi)
K_i <- length(object$coef_names_initial)
object$state_names <- unname(object$state_names)
gamma_pi <- data.frame(
state = unlist(lapply(object$state_names, function(x) x[-1])),
parameter = rep(object$coef_names_initial, each = (S - 1)),
estimate = unlist(eta_pi),
cluster = rep(object$cluster_names, each = (S - 1) * K_i)
state = unlist(object$state_names),
parameter = rep(object$coef_names_initial, each = S),
estimate = unlist(object$gammas$pi),
cluster = rep(object$cluster_names, each = S * K_i)
)
eta_A <- unlist(object$coefficients$eta_A)
K_s <- length(object$coef_names_transition)
gamma_A <- data.frame(
state_from = unlist(object$state_names),
state_to = rep(
unlist(lapply(object$state_names, function(x) x[-1])),
each = S
),
parameter = rep(object$coef_names_transition, each = S * (S - 1)),
estimate = unlist(eta_A),
cluster = rep(object$cluster_names, each = (S - 1) * S * K_s)
state_to = rep(unlist(object$state_names), each = S),
parameter = rep(object$coef_names_transition, each = S * S),
estimate = unlist(object$gammas$A),
cluster = rep(object$cluster_names, each = S * S * K_s)
)
K_o <- length(object$coef_names_emission)
if (object$n_channels == 1) {
gamma_B <- data.frame(
state = unlist(object$state_names),
observations = rep(object$symbol_names[-1], each = S),
parameter = rep(object$coef_names_emission, each = S * (M - 1)),
estimate = unlist(object$coefficients$eta_B),
cluster = rep(object$cluster_names, each = S * (M - 1) * K_o)
observations = rep(object$symbol_names, each = S),
parameter = rep(object$coef_names_emission, each = S * M),
estimate = unlist(object$gammas$B),
cluster = rep(object$cluster_names, each = S * M * K_o)
)
} else {
gamma_B <- data.frame(
state = unlist(object$state_names),
observations = rep(unlist(lapply(object$symbol_names, "[", -1)), each = S),
observations = rep(unlist(object$symbol_names), each = S),
parameter = unlist(lapply(seq_len(object$n_channels), function(i) {
rep(object$coef_names_emission, each = S * (M[i] - 1))
rep(object$coef_names_emission, each = S * M[i])
})),
estimate = unlist(object$coefficients$eta_B),
estimate = unlist(object$gammas$B),
cluster = unlist(lapply(seq_len(object$n_channels), function(i) {
rep(object$cluster_names, each = S * (M[i] - 1) * K_o)
rep(object$cluster_names, each = S * M[i] * K_o)
}))
)
}
eta_omega <- c(object$coefficients$eta_omega)
gamma_omega <- data.frame(
cluster = object$cluster_names[-1],
parameter = rep(object$coef_names_cluster, each = D - 1),
estimate = eta_omega
cluster = object$cluster_names,
parameter = rep(object$coef_names_cluster, each = D),
estimate = c(object$gammas$omega)
)

# stopifnot_(
Expand All @@ -155,10 +145,10 @@ coef.mnhmm <- function(object, probs = c(0.025, 0.5, 0.975), ...) {
# "Argument {.arg probs} must be a {.cls numeric} vector with values
# between 0 and 1."
# )
# p_i <- length(eta_pi)
# p_s <- length(eta_A)
# p_o <- length(eta_B)
# p_d <- length(eta_omega)
# p_i <- length(gamma_pi)
# p_s <- length(gamma_A)
# p_o <- length(gamma_B)
# p_d <- length(gamma_omega)
# sds <- try(
# sqrt(diag(solve(-object$estimation_results$hessian))),
# silent = TRUE
Expand All @@ -173,16 +163,16 @@ coef.mnhmm <- function(object, probs = c(0.025, 0.5, 0.975), ...) {
#
# for(i in seq_along(probs)) {
# q <- qnorm(probs[i])
# gamma_pi[paste0("q", 100 * probs[i])] <- eta_pi + q * sds[seq_len(p_i)]
# gamma_A[paste0("q", 100 * probs[i])] <- eta_A + q * sds[p_i + seq_len(p_s)]
# gamma_B[paste0("q", 100 * probs[i])] <- eta_B + q * sds[p_i + p_s + seq_len(p_o)]
# gamma_omega[paste0("q", 100 * probs[i])] <- eta_omega + q * sds[p_i + p_s + p_o + seq_len(p_d)]
# gamma_pi[paste0("q", 100 * probs[i])] <- gamma_pi + q * sds[seq_len(p_i)]
# gamma_A[paste0("q", 100 * probs[i])] <- gamma_A + q * sds[p_i + seq_len(p_s)]
# gamma_B[paste0("q", 100 * probs[i])] <- gamma_B + q * sds[p_i + p_s + seq_len(p_o)]
# gamma_omega[paste0("q", 100 * probs[i])] <- gamma_omega + q * sds[p_i + p_s + p_o + seq_len(p_d)]
# }

list(
gamma_initial = gamma_pi,
gamma_transition = gamma_A,
gamma_emission = gamma_B,
gamma_cluster = gamma_omega
initial = gamma_pi,
transition = gamma_A,
emission = gamma_B,
cluster = gamma_omega
)
}
6 changes: 3 additions & 3 deletions R/create_base_nhmm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,14 +127,14 @@ create_base_nhmm <- function(observations, data, time, id, n_states,
omega <- model_matrix_cluster_formula(
cluster_formula, data, n_sequences, n_clusters, time, id
)
coefficients <- create_initial_values(
etas <- create_initial_values(
list(gamma_pi = NULL, gamma_A = NULL, gamma_B = NULL, gamma_omega = NULL),
n_states, n_symbols, 0,
length(pi$coef_names), length(A$coef_names), length(B$coef_names),
length(omega$coef_names), n_clusters
)
} else {
coefficients <- create_initial_values(
etas <- create_initial_values(
list(pi = NULL, A = NULL, B = NULL),
n_states, n_symbols, 0,
length(pi$coef_names), length(A$coef_names), length(B$coef_names)
Expand All @@ -154,7 +154,7 @@ create_base_nhmm <- function(observations, data, time, id, n_states,
transition_formula = A$formula,
emission_formula = B$formula,
cluster_formula = if(mixture) omega$formula else NULL,
coefficients = coefficients,
etas = etas,
state_names = state_names,
symbol_names = attr(observations, "symbol_names"),
channel_names = attr(observations, "channel_names"),
Expand Down
Loading

0 comments on commit be42ffa

Please sign in to comment.