R/hmmviterbi2.R

In ziphsmm: Zero-Inflated Poisson Hidden (Semi-)Markov Models

########################################################
#' Viterbi algorithm to decode the latent states in hidden Markov models
#' with covariate values
#' @param y the observed series to be decoded
#' @param ntimes vector specifying the lengths of individual, 
#' i.e. independent, time series. If not specified, the responses are assumed to 
#' form a single time series, i.e. ntimes=length(y)
#' @param M number of latent states
#' @param workparm a vector of values for working parameters, which is the last element returned 
#' from hmmfit() function. This consists the generalized logit of prior probabilities 
#' (except for the 1st state), generalized logit of transition probability matrix
#' (except for the 1st column), the  logit of nonzero structural zero proportions, and 
#' the log poisson means
#' @param zero_init a vector containing structural zero proportions in each state, e.g. set
#' zero_init[i] to be 0 if the i-th state is a regular poisson, and otherwise 1.
#' @param prior_x matrix of covariates for generalized logit of prior probabilites (excluding the 
#' 1st probability). Default to NULL.
#' @param tpm_x matrix of covariates for transition probability matrix (excluding the 1st column).
#' Default to NULL.
#' @param emit_x matrix of covariates for the log poisson means. Default to NULL.
#' @param zeroinfl_x matrix of covariates for the nonzero structural zero proportions. Default to NULL.
#' @param plot whether a plot should be returned
#' @param xlim vector specifying the minimum and maximum on the x-axis in the plot. 
#' Default to NULL.
#' @param ylim vector specifying the minimum and maximum on the y-axis in the plot. 
#' Default to NULL.
#' @param ... further arguments to be passed to the plot() function 
#' @return decoded series of latent states
#' @references Walter Zucchini, Iain L. MacDonald, Roland Langrock. Hidden Markov Models for 
#' Time Series: An Introduction Using R, Second Edition. Chapman & Hall/CRC
#' @examples
#' data(CAT)
#' y <- CAT$activity
#' x <- data.matrix(CAT$night)
#' prior_init <- c(0.5,0.2,0.3)
#' emit_init <- c(10,50,100)
#' zero_init <- c(0.5,0,0)
#' omega <- matrix(c(0.5,0.3,0.2,0.4,0.3,0.3,0.2,0.4,0.4),3,3,byrow=TRUE)
#' fit2 <-  hmmfit(y,rep(1440,3),3,prior_init,omega,
#'      emit_init,zero_init, emit_x=x,zeroinfl_x=x,hessian=FALSE,
#'      method="Nelder-Mead", control=list(maxit=500,trace=1))
#' decode <- hmmviterbi2(y,rep(1440,3),3,fit2$working_parameters,zero_init=c(1,0,0),
#'             emit_x=x,zeroinfl_x=x, plot=TRUE, xlab="time", ylab="count",
#'             xlim=c(0,360),ylim=c(0,200))
#'
#'
#' @useDynLib ziphsmm
#' @importFrom Rcpp evalCpp
#' @export
hmmviterbi2 <- function(y, ntimes=NULL, M, workparm, zero_init,
                        prior_x=NULL, tpm_x=NULL, emit_x=NULL,zeroinfl_x=NULL,
                        plot=FALSE, xlim=NULL, ylim=NULL, ...){
  
  if(floor(M)!=M | M<2) stop("The number of latent states must be an integer greater than or equal to 2!")
  if(length(zero_init)!=M) stop("The dimension of zero_init is not M!")
  if(is.null(ntimes)) ntimes <- length(y)
  state <- rep(NA, sum(ntimes))
  
  tempmat <- matrix(0,nrow=length(y),ncol=1) #for NULL covariates
  if(is.null(prior_x)) {
    ncovprior <- 0
    covprior <- tempmat
  }else{
    ncovprior <- ncol(prior_x)
    covprior <- prior_x
  }
  
  #tpm_x, emit_x, zeroinfl_x
  if(is.null(tpm_x)){
    ncovtpm <- 0
    covtpm <- tempmat
  }else{
    ncovtpm <- ncol(tpm_x)
    covtpm <- tpm_x
  }
  if(is.null(emit_x)){
    ncovemit <- 0
    covemit <- ncol(emit_x)
  }else{
    ncovemit <- ncol(emit_x)
    covemit <- emit_x
  }
  if(is.null(zeroinfl_x)){
    ncovzeroinfl <- 0
    covzeroinfl <- tempmat
  }else{
    ncovzeroinfl <- ncol(zeroinfl_x)
    covzeroinfl <- zeroinfl_x
  }
  

    lastindex <- 0
    for(i in 1:length(ntimes)){
      state[(lastindex+1):(lastindex+ntimes[i])] <- hmm_cov_viterbi(workparm, M,
                                                        y[(lastindex+1):(lastindex+ntimes[i])],
                                                        ncovprior,covprior[(lastindex+1):(lastindex+ntimes[i]),,drop=FALSE],
                                                        ncovtpm,covtpm[(lastindex+1):(lastindex+ntimes[i]),,drop=FALSE],
                                                        ncovzeroinfl,covzeroinfl[(lastindex+1):(lastindex+ntimes[i]),,drop=FALSE],
                                                        ncovemit,covemit[(lastindex+1):(lastindex+ntimes[i]),,drop=FALSE],
                                                        zero_init)
      lastindex <- lastindex+ntimes[i]
    }
  
    if(plot==TRUE){
      xlimit <- rep(NA,2)
      ylimit <- rep(NA,2)
      if(is.null(xlim)){
        xlimit[1] <- 0
        xlimit[2] <- sum(ntimes)
      }else{xlimit <- xlim}
      if(is.null(ylim)){
        ylimit[1] <- 0
        ylimit[2] <- max(y) * 1.3
      }else{ylimit <- ylim}
      
      temp <- y[1:min(length(y),floor(xlimit[2]))]
      plot(temp, xlim=xlimit, ylim=ylimit, type="l",...)
      points(1:length(temp),rep(ylimit[1],length(temp)),
             pch=16,cex=0.8,col=state+1)
      legend <- NULL
      for(i in 1:M) legend <- c(legend,paste("state ",i,sep=""))
      legend("top",legend,pch=16,col=2:(M+1),bty="n",cex=0.9,
             ncol=M,xpd=TRUE)
    }    
    
  return(state)
}