R/mostLikelyHzSequence.R

In aqp: Algorithms for Quantitative Pedology

# what is the most likely sequence, given a markovchain and initial state
# the result isn't likely correct when there are non-zero ties in t
#' @param mc Passed to `markovchain` `conditionalDistribution()`
#' @param t0 Passed to `markovchain` `conditionalDistribution()`
#' @param maxIterations Maximum number of iterations. Default: `10`
#'
#' @export
#' @rdname hzTransitionProbabilities
mostLikelyHzSequence <- function(mc, t0, maxIterations = 10) {

  if(!requireNamespace('markovchain', quietly = TRUE))
    stop('pleast install the `markovchain` package.', call. = FALSE)

  # check for ties
  if(!is.null(attr(as(mc, 'matrix'), 'ties')))
    if(attr(as(mc, 'matrix'), 'ties'))
      warning('ties in transition probability matrix, results may not be reliable', call. = FALSE)

  # store sequence here
  s <- vector(mode = 'character')
  # save the intial state
  i <- 1
  s[i] <- t0
  # compute probabilities for the next state
  cd <- markovchain::conditionalDistribution(mc, t0)

  # search for all but the current state
  not.this.state <- which(!names(cd) == t0)

  ## TODO: this doesn't work when there are ties
  next.state <- names(which.max(cd[not.this.state]))
  i <- i + 1
  s[i] <- next.state
  # continue searching for the most likely next state
  # until reaching the second absorbing state
  # the first absorbing state is retained in 's'
  while(! next.state %in% markovchain::absorbingStates(mc)) {
    cd <- markovchain::conditionalDistribution(mc, next.state)
    # search for all but the current state
    not.this.state <- which(!names(cd) == next.state)
    next.state <- names(which.max(cd[not.this.state]))
    i <- i + 1
    s[i] <- next.state
    # trap run-away iteration
    if(i > maxIterations)
      break
  }
  
  return(s)
}