#' qemiss_gr
#'
#' \code{qemiss_gr} calculates the gradient of \code{qemiss}.
#'
#' @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 vector of length 42, the gradient for qemiss.
#'
#' @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 numeric approximation
#' D1 <- qemiss_gr(pars,X,E,L)
#' require(numDeriv)
#' D2 <- grad(function(u) qemiss(u,X,E,L) , pars)
#' print(round(D1-D2, 10))
qemiss_gr <- function(pars, X, E, L){
pars <- abs(pars)
a <- pars[1:21] # short for alpha
b <- pars[22:42] # short for beta
da <- db <- rep(0, 21)
for (k in 1:21){
for (i in 1:length(X)){
for (t in 1:length(X[[i]])){
da[k] <- da[k] - (sum(1/seq(a[k],a[k]+X[[i]][t]-1)) - log(1+E[i]/b[k]))*L[[i]][t,k]
db[k] <- db[k] - (a[k]/b[k] - (a[k]+X[[i]][t])/(E[i]+b[k])) * L[[i]][t,k]
}
}
}
return(c(da, db))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.