Replacing k_threshold by probs_threshold

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

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
             {