update_ab: update_ab

In shimlab/riboHMM2: Analysis tool for Ribo-seq

Description

update_ab takes an initial parameters c(alpha[k], beta[k]) and update it once.

Usage

update_ab(init, k, X, E, L, ...)

Arguments

init	a vector of length 2, current parameter c(alpha[k], beta[k])
k	a scalar indicating which state is calculated

Value

A vector of length 3, the updated paratemter and fval c(a'[k], b'[k], fval)

Examples

set.seed(2019)
eL <- c(24, 36, 24, 96, 20)
eH <- c(5, 50,35,20, 40,20,10, 50,35,20,
        5, 60,35,10, 55,25,10, 60,35,10, 5)
v <- rep(1000,21); E <- 10; r <- 1

df1 <- generate_uORF(eL, eH, v, E, r)
df2 <- generate_uORF(eL, eH, v, E, r)
la1 <- forwardAlg(df1$x, df1$RNA, df1$trans, df1$v, df1$v/df1$m, df1$E)
lb1 <- backwardAlg(df1$x, df1$RNA, df1$trans, df1$v, df1$v/df1$m, df1$E)
la2 <- forwardAlg(df2$x, df2$RNA, df2$trans, df2$v, df2$v/df2$m, df2$E)
lb2 <- backwardAlg(df2$x, df2$RNA, df2$trans, df2$v, df2$v/df2$m, df2$E)

X <- list(x1=df1$x, x2=df2$x)
L <- list(L1=computeL(la1,lb1), L2=computeL(la2,lb2))
pars <- c(df1$v, df1$v/df1$m)
E <- c(df1$E, df2$E)

k=1
con = list(maxit=1e6, reltol=1e-10)
update_ab(pars[c(k ,21+k)], k, X, E, L, control=con)   # 1000.54292  197.28988   71.11162
update_ab(c(5,1), 1, X, E, L, control=con)             # 793.43812 156.45240  71.12018
update_ab(c(1000 ,5), 1, X, E, L, control=con)         # 5174.59740 1020.34307   71.09436
# 1000.54292 / 197.28988 = 5.071436
# 793.43812 / 156.45240 = 5.071435
# 5174.59740 / 1020.34307 = 5.071429
# optimization of alpha and beta are very unstable, (fval are all close to 71)
# while m = alpha/beta is stable (around 5)

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