diff --git a/R/dd_ML.R b/R/dd_ML.R index 4d59d18..7c1fd4c 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -91,14 +91,15 @@ parsfixdefault = function(ddmodel,brts,missnumspec,idparsopt) #' @param btorph Sets whether the likelihood is for the branching times (0) or #' the phylogeny (1) #' @param soc Sets whether stem or crown age should be used (1 or 2) -#' @param tol Sets the tolerances in the optimization. Consists of: \cr reltolx -#' = relative tolerance of parameter values in optimization \cr reltolf = -#' relative tolerance of function value in optimization \cr abstolx = absolute -#' tolerance of parameter values in optimization -#' @param tolint Sets the tolerance of the numerical integration. COnsists of: -#' \cr absoltint = absolute tolerance and \cr reltolint = relative tolerance. -#' @param k_threshold Sets the threshold of number of species above which -#' logarithmic integration must be used. Default is Inf. +#' @param tol Sets the tolerances in the optimization. Consists of:\cr +#' reltolx = relative tolerance of parameter values in optimization \cr +#' reltolf = relative tolerance of function value in optimization \cr +#' abstolx = absolute tolerance of parameter values in optimization +#' @param tolint Sets the tolerance of the numerical integration. Consists of: \cr +#' absoltint = absolute tolerance and \cr +#' reltolint = relative tolerance. +#' @param probs_threshold Sets the threshold of the probability below which +#' logarithmic integration must be used. Default is 0. #' @param maxiter Sets the maximum number of iterations in the optimization #' @param changeloglikifnoconv if TRUE the loglik will be set to -Inf if ML #' does not converge @@ -150,7 +151,7 @@ dd_ML = function( soc = 2, tol = c(1E-3, 1E-4, 1E-6), tolint = c(1E-10,1E-8), - k_threshold = 100, + probs_threshold = 0, maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = 'simplex', @@ -189,7 +190,7 @@ dd_ML = function( trparsopt[which(initparsopt == Inf)] = 1 trparsfix = parsfix/(1 + parsfix) trparsfix[which(parsfix == Inf)] = 1 - pars2 = c(res,ddmodel,cond,btorph,verbose,soc,tol,maxiter,abstolint = tolint[1],reltolint = tolint[2],k_threshold = k_threshold) + pars2 = c(res,ddmodel,cond,btorph,verbose,soc,tol,maxiter,abstolint = tolint[1],reltolint = tolint[2],probs_threshold = probs_threshold) optimpars = c(tol,maxiter) initloglik = dd_loglik_choosepar(trparsopt = trparsopt,trparsfix = trparsfix,idparsopt = idparsopt,idparsfix = idparsfix,pars2 = pars2,brts = brts,missnumspec = missnumspec, methode = methode) cat("The loglikelihood for the initial parameter values is",initloglik,"\n") diff --git a/R/dd_loglik.R b/R/dd_loglik.R index c1bf988..3fe4318 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -188,7 +188,7 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt } else { if(is.na(pars2['abstolint'])) abstolint <- 1e-10 else abstolint <- pars2['abstolint'] if(is.na(pars2['reltolint'])) reltolint <- 1e-8 else reltolint <- pars2['reltolint'] - if(is.na(pars2['k_threshold'])) k_threshold <- Inf else k_threshold <- pars2['k_threshold'] + if(is.na(pars2['probs_threshold'])) probs_threshold <- 0 else probs_threshold <- pars2['probs_threshold'] brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) if(sum(brts == 0) == 0) { @@ -216,17 +216,17 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt probs = rep(0,lx) probs[1] = 1 # change if other species at stem/crown age lx_old <- lx + logtrafo <- FALSE for(k in 2:(S + 2 - soc)) { k1 = k + (soc - 2) probs <- probs[1:lx] - if(k >= k_threshold) { + if(probs[1] > 0 & probs[1] < probs_threshold) { lx <- min(lx,which(probs == 0) - 1) probs <- probs[1:lx] - if(k == k_threshold) { - rhs_func_name <- 'dd_loglik_log_rhs' - probs <- log(probs) - } + logtrafo <- TRUE + rhs_func_name <- 'dd_loglik_log_rhs' + probs <- log(probs) } y = dd_integrate(probs,brts[(k-1):k],rhs_func_name,c(pars1,k1,ddep),rtol = reltolint,atol = abstolint,method = methode) @@ -239,7 +239,7 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt } if(k < (S + 2 - soc)) { - if(k >= k_threshold) { + if(logtrafo) { fac <- flavec(ddep,la,mu,K,r,lx,k1) lx <- min(lx,which(fac == 0) - 1) probs <- probs[1:lx] @@ -267,12 +267,12 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt cp <- check_probs(loglik,probs[1:lx],verbose); loglik <- cp[[1]]; probs[1:lx] <- cp[[2]]; } } - if((k < k_threshold & probs[1 + missnumspec] <= 0) | loglik == -Inf | is.na(loglik) | is.nan(loglik)) + if((!logtrafo & probs[1 + missnumspec] <= 0) | loglik == -Inf | is.na(loglik) | is.nan(loglik)) { if(verbose) cat('Probabilities smaller than 0 or other numerical problems are encountered in final result.\n') loglik = -Inf } else { - if(k < k_threshold) { + if(!logtrafo) { loglik = loglik + (cond != 3 | soc == 1) * log(probs[1 + (cond != 3) * missnumspec]) - lgamma(S + missnumspec + 1) + lgamma(S + 1) + lgamma(missnumspec + 1) } else { diff --git a/man/dd_ML.Rd b/man/dd_ML.Rd index 3ed2c67..198514a 100644 --- a/man/dd_ML.Rd +++ b/man/dd_ML.Rd @@ -19,7 +19,7 @@ dd_ML( soc = 2, tol = c(0.001, 1e-04, 1e-06), tolint = c(1e-10, 1e-08), - k_threshold = 100, + probs_threshold = 0, maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "simplex", @@ -97,16 +97,17 @@ the phylogeny (1)} \item{soc}{Sets whether stem or crown age should be used (1 or 2)} -\item{tol}{Sets the tolerances in the optimization. Consists of: \cr reltolx -= relative tolerance of parameter values in optimization \cr reltolf = -relative tolerance of function value in optimization \cr abstolx = absolute -tolerance of parameter values in optimization} +\item{tol}{Sets the tolerances in the optimization. Consists of:\cr +reltolx = relative tolerance of parameter values in optimization \cr +reltolf = relative tolerance of function value in optimization \cr +abstolx = absolute tolerance of parameter values in optimization} -\item{tolint}{Sets the tolerance of the numerical integration. COnsists of: -\cr absoltint = absolute tolerance and \cr reltolint = relative tolerance.} +\item{tolint}{Sets the tolerance of the numerical integration. Consists of: \cr +absoltint = absolute tolerance and \cr +reltolint = relative tolerance.} -\item{k_threshold}{Sets the threshold of number of species above which -logarithmic integration must be used. Default is Inf.} +\item{probs_threshold}{Sets the threshold of the probability below which +logarithmic integration must be used. Default is 0.} \item{maxiter}{Sets the maximum number of iterations in the optimization}