R/qemiss.R

In shimlab/riboHMM2: Analysis tool for Ribo-seq

#' qemiss
#'
#' \code{qemiss} is the emission part of the Q function in E-step,
#' i.e.: Sum_{i,t,k} log[fk(x_t)] * L^i_k(t). (with negative sign)
#'
#' @export
#' @param pars a vector of length 42. c(alpha, beta)
#' @param X    a list of vectors of observed states x
#' @param E    a vector of normalizing constant for each observed chain in X
#' @param L    a list of matrix L from \code{computeL}
#' @return     A scalar, the (negative) value of the target function
#'                that would later be minimized.
#'
#' @examples
#' df <- uORF
#' X <- L <- list()
#' E <- c()
#' for (i in 1:2){
#'   X[[i]] <- df[[i]]$x
#'   RNA <- df[[i]]$RNA
#'   E[i]=df[[i]]$E;   trans=df[[i]]$trans;
#'   a=df[[i]]$v;      b=df[[i]]$v/df[[i]]$m
#'   la <- forwardAlg(X[[i]], RNA, trans, a, b, E[i])
#'   lb <- backwardAlg(X[[i]], RNA, trans, a, b, E[i])
#'   L[[i]] <- computeL(la, lb)
#' }
#' pars <- c(df[[1]]$v, df[[1]]$v/df[[1]]$m)
#'
#' qemiss(pars,X,E,L)

qemiss <- function(pars, X, E, L){
  pars <- abs(pars)
  alpha <- pars[1:21]
  beta <- pars[22:42]
  qe = 0
  for (i in 1:length(X)){
    for (t in 1:length(X[[i]])){
      lf_t <- lnb(X[[i]][t], alpha, beta, E[i])
      qe <- qe + sum(lf_t * L[[i]][t,])
    }
  }
  return(-qe)
}