Nothing
#' Posterior decoding.
#'
#' Calculate the posterior probability of a sequence given a model.
#'
#' @param x an object of class \code{'HMM'} or \code{'PHMM'}.
#' @param y a vector of mode "character" or "raw" (a "DNAbin" or "AAbin"
#' object) representing a single sequence hypothetically emitted by
#' the model in \code{x}.
#' @param logspace logical indicating whether the emission and transition
#' probabilities of x are logged. If \code{logspace = "autodetect"}
#' (default setting), the function will automatically detect
#' if the probabilities are logged, returning an error if
#' inconsistencies are found. Note that choosing the latter option
#' increases the computational overhead; therefore specifying
#' \code{TRUE} or \code{FALSE} can reduce the running time.
#' @param cpp logical, indicates whether the dynamic programming matrix
#' should be filled using compiled C++ functions (default; many times faster).
#' The R version is primarily retained for bug-fixing and experimentation.
#' @param ... additional arguments to be passed between methods.
#' @return a vector, matrix or array of posterior probabilities.
#' @details
#' See Durbin et al (1998) chapter 3.2 for details on the calculation
#' and interpretation of posterior state probabilities.
#' Currently no method is available for profile HMMs, but this may be
#' included in a future version if required.
#' @author Shaun Wilkinson
#' @references
#' Durbin R, Eddy SR, Krogh A, Mitchison G (1998) Biological
#' sequence analysis: probabilistic models of proteins and nucleic acids.
#' Cambridge University Press, Cambridge, United Kingdom.
#' @examples
#' ## Posterior decoding for standard hidden Markov models
#' ## The dishonest casino example from Durbin et al (1998) chapter 3.2
#' states <- c("Begin", "Fair", "Loaded")
#' residues <- paste(1:6)
#' ### Define the transition probability matrix
#' A <- matrix(c(0, 0, 0, 0.99, 0.95, 0.1, 0.01, 0.05, 0.9), nrow = 3)
#' dimnames(A) <- list(from = states, to = states)
#' ### Define the emission probability matrix
#' E <- matrix(c(rep(1/6, 6), rep(1/10, 5), 1/2), nrow = 2, byrow = TRUE)
#' dimnames(E) <- list(states = states[-1], residues = residues)
#' ### Build and plot the HMM object
#' x <- structure(list(A = A, E = E), class = "HMM")
#' plot(x, main = "Dishonest casino HMM")
#' ### Calculate posterior probabilities
#' data(casino)
#' casino.post <- posterior(x, casino)
#' plot(1:300, casino.post[1, ], type = "l", xlab = "Roll number",
#' ylab = "Posterior probability of dice being fair",
#' main = "The dishonest casino")
#' @seealso \code{\link{forward}}, \code{\link{backward}}, \code{\link{Viterbi}}
#' @name posterior
################################################################################
posterior <- function(x, y, ...){
UseMethod("posterior")
}
################################################################################
#' @rdname posterior
################################################################################
posterior.HMM <- function(x, y, logspace = "autodetect", cpp = TRUE, ...){
if(identical(logspace, 'autodetect')) logspace <- .logdetect(x)
back <- backward(x, y, logspace = logspace, cpp = cpp)
B <- back$array
forw <- forward(x, y, logspace = logspace, cpp = cpp)
R <- forw$array
logPx <- forw$score
postprobs <- exp(R + B - logPx)
return(postprobs)
}
################################################################################
#' @rdname posterior
################################################################################
posterior.PHMM <- function(x, y, logspace = "autodetect", cpp = TRUE, ...){
### stop("Posterior method not currently available for profile HMMs")
### placeholder
if(identical(logspace, 'autodetect')) logspace <- .logdetect(x)
back <- backward(x, y, logspace = logspace, cpp = cpp)
B <- back$array
forw <- forward(x, y, logspace = logspace, cpp = cpp)
R <- forw$array
logPx <- forw$score
postprobs <- exp(R + B - logPx)
return(postprobs)
}
################################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.