R/qemiss_gr2.R

In shimlab/riboHMM2: Analysis tool for Ribo-seq

#' qemiss_gr2
#'
#' \code{qemiss_gr2} calculates the gradient of \code{qemiss2}.
#' \code{qemiss_gr2} calculates the target function for each k
#' instead of all k.
#'
#' @export
#' @param pars a vector of length 2. c(alpha[k], beta[k])
#' @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}
#' @param k    a scalar indicating which state is calculated
#' @return     A vector of length 2, the gradient for qemiss2.
#'
#' @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)
#'
#' # check by comparing with qemiss_gr & numeric approximation
#' D1 <- qemiss_gr(pars,X,E,L)
#' require(numDeriv)
#' D2 <- grad(function(u) qemiss(u,X,E,L) , pars)
#' D3 <- rep(0,42)
#' for (k in 1:21){
#'   D3[c(k,21+k)] <- qemiss_gr2(pars[c(k,21+k)], X, E, L, k)
#' }
#' print(round(D1-D2, 10))
#' print(round(D3-D1, 10))
#' print(round(D3-D2, 10))

qemiss_gr2 <- function(pars, X, E, L, k){
  pars <- abs(pars)
  a <- pars[1]    # short for alpha
  b <- pars[2]   # short for beta
  dak <- dbk <- 0
  for (i in 1:length(X)){
    for (t in 1:length(X[[i]])){
      if (X[[i]][t]==0){
        dak <- dak - (-log(1+E[i]/b))*L[[i]][t,k]
      }else{
        s <- sum(1/seq(a,a+X[[i]][t]-1))
        dak <- dak - (s - log(1+E[i]/b))*L[[i]][t,k]
      }
      dbk <- dbk - (a/b - (a+X[[i]][t])/(E[i]+b)) * L[[i]][t,k]
    }
  }
  return(c(dak, dbk))
}