R/computeH.R

In shimlab/riboHMM2: Analysis tool for Ribo-seq

#' computeH
#'
#' \code{computeH} calculates H_kj (t) = P( Z[t]=k, Z[t+1]=j | x[1:n] )
#'
#' @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.
#' @param la    a matrix from Forward Algorithm
#' @param lb    a matrix from Backward Algorithm
#' @return      A matrix H
#'
#' @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)
#' lb <- backwardAlg(x, RNA, trans, a, b, E)
#' L <- computeL(la, lb)
#' H <- computeH(x, RNA, trans, a, b, E, la, lb)
#'
#' H <- round(H, 2)
#' View(H)
#' df$z

computeH <- function(x, RNA, trans, alpha, beta, E, la, lb){
  n <- length(x)
  la_n <- la[n,]
  lse <- logSumExp(la_n)
  n_pars <- 5
  H <- matrix(0, n-1, n_pars)
  kk <- c(1,1,1,11,11)
  jj <- c(1,2,12,11,12)
  for (t in 1:(n-1)){
    lA <- log(transprob(trans, RNA[t], RNA[t]))
    for (i in 1:n_pars){
      lf <- lnb(x[t+1], alpha[jj[i]], beta[jj[i]], E)
      H[t,i] <- exp(la[t,kk[i]] + lb[t+1,jj[i]] + lA[kk[i],jj[i]] + lf - lse)
    }
  }
  colnames(H) <- c('1,1', '1,2', '1,12', '11,11','11,12')
  return(H)
}