R/forwardAlg.R

In shimlab/riboHMM2: Analysis tool for Ribo-seq

#' forwardAlg
#'
#' \code{forwardAlg} performs the forward algorithm on a single trajecture of observed chain x.
#' \bold{a} (in log scale), where \bold{a[k,t] = P(Z[t]=k, x[1:t])}
#'
#' @export
#' @param x     a vector of observed states
#' @param RNA   a 0-1 vector. 1 if next 3-base is stop codon
#' @param trans a vector c(rho_u, rho, delta)
#' @param alpha shape parameter in gamma distribution
#' @param beta  rate parameter in gamma distribution
#' @param E     a scalar. Normalizing constant for the observed chain x.
#' @return      A matrix, \bold{a} (in log scale).
#'
#' @examples
#' df <- uORF[[1]]
#' x=df$x; RNA=df$RNA; trans=df$trans; a=df$v; b=df$v/df$m; E=df$E
#'
#' la <- forwardAlg(x, RNA, trans, a, b, E)

forwardAlg <- function(x, RNA, trans, alpha, beta, E){
  n <- length(x)
  M <- 21
  la <- matrix(-Inf, n, M)
  la[1,1] <- lnb(x[1], alpha[1], beta[1], E)

  for (t in 2:n){
    lA <- log(transprob(trans, RNA[t-1], RNA[t-1]))
    for (k in 1:M){
      ls <- la[t-1,] + lnb(x[t], alpha[k], beta[k], E) + lA[,k]
      la[t,k] <- logSumExp(ls)
    }
  }
  colnames(la) <- seq(1,21)
  return(la)
}