#' q_mv_gr2
#'
#' \code{q_mv_gr2} calculates the gradient of \code{q_mv2}.
#' \code{q_mv_gr2} calculates the target function for each k
#' instead of all k.
#'
#' @export
#' @param pars a vector of length 2. c(v[k], m[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 q_mv2.
#'
#' @examples
#' df <- uORF[1:10]
#' X <- L <- list()
#' E <- c()
#' for (i in 1:length(df)){
#' X[[i]] <- df[[i]]$x
#' RNA[[i]] <- 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[[i]], trans, a, b, E[i])
#' lb <- backwardAlg(X[[i]], RNA[[i]], trans, a, b, E[i])
#' L[[i]] <- computeL(la, lb)
#' }
#' pars <- c(df[[1]]$v, df[[1]]$m)
#'
#' # check by comparing with numeric approximation
#' require(numDeriv)
#' D1 <- rep(0,42)
#' for (k in 1:21){
#' D1[c(k,21+k)] <- grad(function(u) q_mv2(u,X,E,L,k) , pars[c(k,21+k)])
#' }
#' D2 <- rep(0,42)
#' for (k in 1:21){
#' D2[c(k,21+k)] <- q_mv_gr2(pars[c(k,21+k)], X, E, L, k)
#' }
#' print(round(D1-D2, 10))
q_mv_gr2 <- function(pars, X, E, L, k){
pars <- abs(pars)
v <- pars[1]
m <- pars[2]
dvk <- dmk <- 0
for (i in 1:length(X)){
for (t in 1:length(X[[i]])){
if (X[[i]][t]==0){
dvk <- dvk - (-log(1+m*E[i]/v) + (m*E[i]-X[[i]][t])/(m*E[i]+v)) * L[[i]][t,k]
}else{
s <- sum(1/seq(v,v+X[[i]][t]-1))
dvk <- dvk - (s - log(1+m*E[i]/v) + (m*E[i]-X[[i]][t])/(m*E[i]+v)) * L[[i]][t,k]
}
dmk <- dmk - (v/m * (X[[i]][t]-m*E[i])/(m*E[i]+v)) * L[[i]][t,k]
}
}
return(c(dvk, dmk))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.