#' Forward and Backward Probabilities for Hidden Markov Model
#'
#' The \code{forward_backward} function computes scaled forward and backward probabilities of a hidden Markov model.
#'
#' @export
#' @param model Object of class \code{hmm} or \code{mhmm}.
#' @param forward_only If \code{TRUE}, only forward probabilities are computed. The default is \code{FALSE}.
#' @param log_space Compute forward and backward probabilities in logarithmic scale instead of scaling.
#' The default is \code{FALSE}.
#' @param threads Number of threads used in parallel computing. The default is 1.
#' @return List with components
#' \item{forward_probs}{If \code{log_space = FALSE}, scaled forward probabilities, i.e. probability of state given
#' observations up to that time point. If \code{log_space = TRUE},
#' logarithms of non-scaled forward probabilities. }
#' \item{backward_probs}{Scaled backward probabilities (\code{log_space = FALSE}),
#' or logarithms of non-scaled backward probabilities(\code{log_space = TRUE}). }
#' \item{scaling_factors}{Sum of non-scaled forward probabilities at each time point.
#' Only computed if \code{log_space = FALSE}.}
#' In case of multiple observations, these are computed independently for each sequence.
#' @examples
#' # Load a pre-defined MHMM
#' data("mhmm_biofam")
#'
#' # Compute forward and backward probabilities
#' fb <- forward_backward(mhmm_biofam)
#'
#' # The most probable hidden state at time t
#' # given the observations up to time t for the first subject:
#' apply(fb$forward_probs[, , 1], 2, which.max)
#'
forward_backward <- function(model, forward_only = FALSE, log_space = FALSE, threads = 1) {
if (!inherits(model, c("hmm", "mhmm"))) {
stop("Argument model must be an object of class 'hmm' or 'mhmm.")
}
if (inherits(model, "mhmm")) {
mix <- TRUE
model <- combine_models(model)
} else {
mix <- FALSE
}
if (model$n_channels == 1) {
model$observations <- list(model$observations)
model$emission_probs <- list(model$emission_probs)
}
obsArray <- array(0, c(model$n_sequences, model$length_of_sequences, model$n_channels))
for (i in 1:model$n_channels) {
obsArray[, , i] <- sapply(model$observations[[i]], as.integer) - 1L
obsArray[, , i][obsArray[, , i] > model$n_symbols[i]] <- model$n_symbols[i]
}
obsArray <- aperm(obsArray)
emissionArray <- array(1, c(model$n_states, max(model$n_symbols) + 1, model$n_channels))
for (i in 1:model$n_channels) {
emissionArray[, 1:model$n_symbols[i], i] <- model$emission_probs[[i]]
}
if (mix) {
if (!log_space) {
out <- forwardbackwardx(
model$transition_probs, emissionArray,
model$initial_probs, obsArray, model$coefficients, model$X, model$n_states_in_clusters, forward_only, threads
)
} else {
out <- log_forwardbackwardx(
model$transition_probs, emissionArray,
model$initial_probs, obsArray, model$coefficients, model$X, model$n_states_in_clusters, forward_only, threads
)
}
} else {
if (!log_space) {
out <- forwardbackward(
model$transition_probs, emissionArray,
model$initial_probs, obsArray, forward_only, threads
)
} else {
out <- log_forwardbackward(
model$transition_probs, emissionArray,
model$initial_probs, obsArray, forward_only, threads
)
}
}
if (is.null(time_names <- colnames(model$observations[[1]]))) {
time_names <- 1:model$length_of_sequences
}
if (is.null(sequence_names <- rownames(model$observations[[1]]))) {
sequence_names <- 1:model$n_sequences
}
dimnames(out$forward_probs) <- list(
"state" = model$state_names,
"time" = time_names, "sequence" = sequence_names
)
if (!log_space) {
dimnames(out$scaling_factors) <- list("time" = time_names, "sequence" = sequence_names)
}
if (!forward_only) {
dimnames(out$backward_probs) <- list(
"state" = model$state_names,
"time" = time_names, "sequence" = sequence_names
)
}
out
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.