R/update_mv.R

In shimlab/riboHMM2: Analysis tool for Ribo-seq

#' update_mv
#'
#' \code{update_mv} takes an initial parameters c(v[k], m[k]) and update it once.
#'
#' @export
#' @param init a vector of length 2, current parameter c(v[k], m[k])
#' @param k    a scalar indicating which state is calculated
#' @return     A vector of length 3, the updated paratemter and fval c(v'[k], m'[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))
#' E <- c(df1$E, df2$E)
#'
#' pars <- c(df1$v, df1$m)
#' # con = list(maxit=1e6, reltol=1e-10)  take long time to run
#' for (k in 1:21){
#'   update_mv(pars[c(k ,21+k)], k, X, E, L)
#' }
#'
#' # optimization of v is still unstable since fval does not change much along v
#' k=2
#' update_mv(pars[c(k ,21+k)], k, X, E, L)
#' update_mv(c(500,50), k, X, E, L)
#' update_mv(c(1500 ,50), k, X, E, L)
#'
#' # modify argument \code{control} to get more accurate result
#' con = list(maxit=1e4, reltol=1e-9)
#' update_mv(pars[c(k ,21+k)], k, X, E, L, control= con)
#' update_mv(c(500,5), k, X, E, L, control= con)
#' update_mv(c(1500,5), k, X, E, L, control= con)

update_mv <- function(init, k, X, E, L, ...){
  opt <- optim(init, fn = q_mv2, gr = q_mv_gr2,
               X=X, E=E, L=L, k=k, method="BFGS", ...)
  vm <- opt$par
  fval <- opt$value
  return(c(abs(vm), fval))
}