qemiss_gr2: qemiss_gr2

In shimlab/riboHMM2: Analysis tool for Ribo-seq

Description

qemiss_gr2 calculates the gradient of qemiss2. qemiss_gr2 calculates the target function for each k instead of all k.

Usage

qemiss_gr2(pars, X, E, L, k)

Arguments

pars	a vector of length 2. c(alpha[k], beta[k])
X	a list of vectors of observed states x
E	a vector of normalizing constant for each observed chain in X
L	a list of matrix L from computeL
k	a scalar indicating which state is calculated

Value

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))

shimlab/riboHMM2 documentation built on May 19, 2019, 6:23 p.m.