R/optimalPath.R

In shimlab/riboHMM2: Analysis tool for Ribo-seq

#' 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     vector of single observed trajectures
#' @param RNA   0-1 vector. 1 if next 3-base is stop codon
#' @param E     scalar. Normalizing constant for the observed chain x.
#' @param v     vector of shape parameter in gamma distribution (from \code{estimate()})
#' @param m     vector of mean parameter in gamma distribution (from \code{estimate()})
#' @param trans vector c(rho_u, rho, delta) (from \code{estimate()})
#' @return      List containing:
#'                \itemize{
#'                  \item the optimal path \code{optPth}
#'                  \item the corresponding log-likelihood \code{loglik}
#'                }
#'
#' @examples
#' k <- 5  # number of transcripts from each group (see data ?uORF)
#' ii <- c()
#' for (i in 1:5){
#'   ii <- c(ii, ((i-1)*100+1):((i-1)*100+k))
#' }
#' print(ii)
#'
#' # estimate parameters
#' df <- uORF[ii]
#' num <- length(df)
#' X <- RNA <- list(); E <- c()
#' for (i in 1:num){
#'   X[[i]] <-  df[[i]]$x
#'   RNA[[i]] <- df[[i]]$RNA
#'   E[i] <- df[[i]]$E
#' }
#' init <- init_generate(df, cutoff=10, r=4)
#' est <- estimate(X, RNA, E, init, tol=0.01)
#' print(est)
#'
#' # using estimated parameters in \code{optimalPath}
#' res <- matrix(0, num ,3)
#' colnames(res) <- c("rate", "TRUE", "FASLE")
#' for (i in 1:num){
#'   n <- length(df[[i]]$x)
#'   optPth <- optimalPath(df[[i]]$x, df[[i]]$RNA, E[i], est$v_hat, est$m_hat, est$trans_hat)$optPth
#'   res[i,2] <- sum(optPth == df[[i]]$z)
#'   res[i,3] <- n - res[i,2]
#'   res[i,1] <- res[i,2]/n
#' }
#' print(res)
#'

optimalPath <- function(x, RNA, E, v, m, trans){
  M <- 21
  n <- length(x)
  G <- rep(0, M)
  S <- matrix(0, n-1,M)
  beta <- v/m

  G <- c(0, rep(-Inf, M-1))
  for (t in 2:n){
    lA <- log(transprob(trans, RNA[t-1], RNA[t-1]))
    lf <- lnb(x[t], v, beta, E)
    G_temp <- rep(0, M)
    for (k in 1:M){
      Gk_temp <- G + lA[,k] + lf[k]
      #for (j in 1:M)  Gk_temp[j] <- G[j] + lA[j,k] + lf[k]
      G_temp[k] <- max(Gk_temp)
      S[t-1,k] <- which.max(Gk_temp)
    }
    G <- G_temp
  }

  s <- which.max(G)
  for (t in n:2){
    s <- c(S[t-1, s[1]], s)
  }

  return(list(optPth=s, loglik=max(G)))
}