R/optimalPath.R

In jiangrongo/HMM:

#' optimalPath
#'
#' \code{optimalPath} performs Viterbi Algorithm on a single
#' observed trajectures \code{x}
#' to predict the optimal path
#' according to the rule, Maximum A Posteriori (MAP).
#'
#' @export
#' @param x     a single observed trajectures
#' @param u     (All parameters used are from the result of \code{estimateHMM}) a vector of estimated means of
#'                 Normal distribution for each state (emission probability)
#' @param sig   a vector of estimated standard
#'                 deviations of Normal distribution for each state (emission probability)
#' @param A     a matrix of estimated transition probability
#' @param pi0   a vector of estimated initial probability.
#' @return      A list containing:
#'                \itemize{
#'                  \item the optimal path \code{opt.path}
#'                  \item the corresponding log-likelihood \code{loglik}
#'                }
#'
#' @examples
#' set.seed(1221)
#' num = 10; n = 250
#' df <- generateHMM(num=num,n=n)
#'
#' ### using Constrain 1
#' temp <- estimateHMM(df$X, M=3, constr=matrix(1,3,3), tol=0.001)
#' u <- temp$u[,temp$iter+1]; sig <- temp$sig[,temp$iter+1]
#' A <- temp$trans[[temp$iter+1]]; pi0 <- temp$pi0
#' res <- matrix(0, num ,3)
#' colnames(res) <- c("rate", "TRUE", "FASLE")
#' for (i in 1:num){
#'   res[i,2] <- sum(optimalPath(df$X[,i], u, sig, A, pi0)$opt.path == df$Z[,i])
#'   res[i,3] <- n-res[i,2]
#'   res[i,1] <- res[i,2]/n
#' }
#' print(res)
#'
#'
#' ### using Constrain 2
#' constr2 <- matrix(1,3,3)
#' constr2[1,1] <- 0; constr2[3,3] <- 0;
#' temp <- estimateHMM(df$X, M=3, constr=constr2, tol=0.001)
#' u <- temp$u[,temp$iter+1]; sig <- temp$sig[,temp$iter+1]
#' A <- temp$trans[[temp$iter+1]]; pi0 <- temp$pi0
#' res <- matrix(0, num ,3)
#' colnames(res) <- c("rate", "TRUE", "FASLE")
#' for (i in 1:num){
#'   res[i,2] <- sum(optimalPath(df$X[,i], u, sig, A, pi0)$opt.path == df$Z[,i])
#'   res[i,3] <- n-res[i,2]
#'   res[i,1] <- res[i,2]/n
#' }
#' print(res)

optimalPath <- function(x, u, sig, A, pi0){
  logA <- log(A)
  logpi <- log(pi0)
  M <- length(u)
  tn <- length(x)
  opt.temp <- matrix(1, M, tn-1)
  Gbf <- rep(0, M)
  for (cr in 1:M){
    Gbf[cr] <- logpi[cr] + dnorm(x[1], mean=u[cr], sd=sig[cr], log=TRUE)
  }
  for (t in 2:tn){
    G <- rep(0, M)
    for (cr in 1:M){
      G[cr] <- Gbf[1] + logA[1,cr] + dnorm(x[t], mean=u[cr], sd=sig[cr], log=TRUE)
      for (bf in 2:M){
        temp <- Gbf[bf] + logA[bf,cr] + dnorm(x[t], mean=u[cr], sd=sig[cr], log=TRUE)
        if (temp > G[cr]){
          G[cr] <- temp
          opt.temp[cr, t-1] <- bf
        }
      }
    }
    Gbf <- G
  }
  opt <- c(which.max(G))
  for (t in (tn-1):1){
    opt <- c(opt.temp[head(opt,1),t] ,opt)
  }
  return(list(opt.path=opt, loglik=max(G)))
}