R/controlEM.R

In plbaldoni/mixNBHMM: Differential Peak Calling via Hidden Markov Model With Mixture of Negative Binomial Distributions

#' Control parameters of the EM algorithm from mixNBHMM
#'
#' This function passes controlling parameters for mixNBHMM. Most of these parameters control the EM algorithm.
#'
#' @param epsilon.em Either a positive value or a vector of size 4 with the convergence tolerance values for the EM algorithm (see 'criterion' below). Default is c(1e-3,1e-3,1e-6,1e-3)
#' @param maxit.em integer giving the maximum number of EM iterations (default 500)
#' @param minit.em integer giving the minimum number of EM iterations to start evaluating the convergence (default 3)
#' @param gap.em integer giving the number of EM iterations apart to compute the convergence criterion (default 3)
#' @param maxcount.em integer giving the number of consecutive EM iterations satisfying the convergence criterion  in order to stop the algorithm (default 3)
#' @param max.phi maximum positive value allowed for the dispersion parameters (default 1000)
#' @param criterion convergence criterion: either "MRCPE" (maximum absolute relative change in parameter estimates), "MACPE" (maximum absolute change of parameter estimates),
#' "ARCEL" (absolute relative change of the Q-function), "ACC" (agreement of Viterbi peak calls),
#' or "MULTI" (simultaneously check for MRCPE, MACPE, ARCEL, and ACC).
#' For ACC, it computes the percentage of windows falling in the main diagonal of a 3 by 3 table of Viterbi predictions 'gap.em' iterations apart. Default is "MULTI"
#' @param min.zero minimum positive value allowed in computations to avoid having zeros (default is .Machine$double.xmin)
#' @param pcut cutoff for rejection controlled EM algorithm (default 0.05)
#' @param quiet whether to print messages (default F)
#' @param maxit.innerem integer giving the maximum number of inner EM iterations (default 5)
#' @param epsilon.innerem a positive value with the convergence tolerance value for the inner EM algorithm. The criterion for the inner EM is "MRCPE".
#' @param trim.offset either NULL or an integer indicating the number of decimal places to be used in the offset (default 3)
#' @param pattern either 'all', 'cluster', or a list with length equal to the number of differential patterns to be modeled by the differential HMM state. If 'all',
#' every possible combinatorial pattern will be considered. If 'cluster', the algorithm will determine the number of mixture components based on the number of most representative combinatorial patterns while discarding patterns whose prevalence is less than 0.05. If a list, elements of it should specify the differential patterns to be modeled by each mixture component.
#' For instance, if pattern = list(2,c(1,3)) the mixture model will have two components that will represent the enrichment of condition 2 alone and the enrichment of conditions 1 and 3 together.
#'
#' @return A list with components equal to the arguments
#'
#' @author Pedro L. Baldoni, \email{pedrobaldoni@gmail.com}
#' @references \url{https://github.com/plbaldoni/mixNBHMM}
#'
#' @examples
#' # Setting maxit.em = 100 (no more than 100 EM iterations)
#' control = controlEM(maxit.em = 100)
#'
#' @useDynLib mixNBHMM
#' @importFrom Rcpp evalCpp
#' @export
#'
controlEM = function(epsilon.em=c(1e-4,1e-4,1e-6,1e-4),maxit.em=500,minit.em=3,gap.em=3,maxcount.em=3,max.phi=1e3,criterion='MULTI',
                        min.zero=.Machine$double.xmin,pcut=0.05,quiet=T,maxit.innerem = 5,epsilon.innerem = 1e-3,
                        trim.offset = 3,pattern = 'all'){
    #if (!is.numeric(epsilon.em) || epsilon.em <= 0){stop("value of 'epsilon.em' must be > 0")}
    if (!maxit.em%%1==0 || maxit.em <= 0){stop("value of 'maxit.em' must be a positive integer")}
    if (!minit.em%%1==0 || minit.em <= 0){stop("value of 'minit.em' must be a positive integer")}
    if (!gap.em%%1==0 || gap.em <= 0 || gap.em>minit.em){stop("value of 'gap.em' must be a positive integer <= minit.em")}
    if (!maxcount.em%%1==0 || maxcount.em <= 0){stop("value of 'maxcount.em' must be a positive integer")}
    if (!is.numeric(max.phi) || max.phi <= 0){stop("value of 'max.phi' must be > 0")}
    list(epsilon.em=epsilon.em,min.zero=min.zero,maxit.em=maxit.em,minit.em=minit.em,gap.em=gap.em,
         maxcount.em=maxcount.em,max.phi=max.phi,criterion=criterion,pcut=pcut,quiet=quiet,
         maxit.innerem=maxit.innerem,epsilon.innerem=epsilon.innerem,trim.offset=trim.offset,pattern = pattern)
}