Nothing
#' Loglikelihood for diversity-dependent diversification models with decoupling
#' of a subclade from a main clade at time t = t_d
#'
#' This function computes loglikelihood of a diversity-dependent
#' diversification model for a given set of branching times and parameter
#' values where the diversity-dependent dynamics of a subclade decouple from
#' the dynamics of the main clade at time t_d, potentially accompanied by a
#' shift in parameters.
#'
#'
#' @param pars1 Vector of parameters: \cr \cr
#' \code{pars1[1]} corresponds to
#' lambda_M (speciation rate) of the main clade \cr
#' \code{pars1[2]} corresponds
#' to mu_M (extinction rate) of the main clade \cr
#' \code{pars1[3]} corresponds
#' to K_M (clade-level carrying capacity) of the main clade \cr
#' \code{pars1[4]} corresponds to lambda_S (speciation rate) of the subclade \cr
#' \code{pars1[5]} corresponds to mu_S (extinction rate) of the subclade \cr
#' \code{pars1[6]} corresponds to K_S (clade-level carrying capacity) of the
#' subclade \cr
#' \code{pars1[7]} corresponds to t_d (the time of decoupling)
#' @param pars2 Vector of model settings: \cr \cr
#' \code{pars2[1]} sets the
#' maximum number of species for which a probability must be computed. This
#' must be larger than 1 + missnumspec + length(brts). \cr \cr
#' \code{pars2[2]} sets the model of diversity-dependence: \cr
#' - \code{pars2[2] == 1} linear
#' dependence in speciation rate with parameter K (= diversity where speciation
#' = extinction)\cr
#' - \code{pars2[2] == 1.3} linear dependence in speciation
#' rate with parameter K' (= diversity where speciation = 0)\cr
#' - \code{pars2[2] == 2} exponential dependence in speciation rate with
#' parameter K (= diversity where speciation = extinction)\cr
#' - \code{pars2[2] == 2.1} variant of exponential dependence in speciation rate with offset at
#' infinity\cr
#' - \code{pars2[2] == 2.2} 1/n dependence in speciation rate\cr
#' - \code{pars2[2] == 2.3} exponential dependence in speciation rate with
#' parameter x (= exponent)\cr
#' - \code{pars2[2] == 3} linear dependence in
#' extinction rate \cr
#' - \code{pars2[2] == 4} exponential dependence in
#' extinction rate \cr
#' - \code{pars2[2] == 4.1} variant of exponential
#' dependence in extinction rate with offset at infinity\cr
#' - \code{pars2[2] == 4.2} 1/n dependence in extinction rate\cr\cr
#' \code{pars2[3]} sets the
#' conditioning: \cr
#' - \code{pars2[3] == 0} no conditioning (or just crown age) \cr
#' - \code{pars2[3] == 1} conditioning on non-extinction of the phylogeny \cr
#' - \code{pars2[3] == 2} conditioning on number of species and crown age;
#' not yet implemented \cr
#' - \code{pars2[3] == 3} conditioning on number of species only;
#' not yet implemented \cr
#' - \code{pars2[3] == 4} conditioning on survival of the subclade \cr
#' - \code{pars2[3] == 5} conditioning on survival of all subclades
#' and of both crown lineages in the main clade. This assumes that subclades
#' that have already shifted do not undergo another shift, i.e. shifts only
#' occur in the main clade. \cr \cr
#' \code{pars2[4]} Obsolete. \cr \cr
#' \code{pars2[5]} sets whether the parameters and likelihood should be
#' shown on screen (1) or not (0) \cr \cr
#' \code{pars2[6]} sets whether the first data point is stem age (1) or crown age (2) \cr\cr
#' \code{pars2[7]} sets whether the old (incorrect) likelihood should be used (0),
#' or whether the new corrected likelihood should be used (1).
#' @param brtsM A set of branching times of the main clade in the phylogeny,
#' all positive
#' @param brtsS A set of branching times of the subclade in the phylogeny, all
#' positive
#' @param missnumspec The number of species that are in the clade but missing
#' in the phylogeny. One can specify the sum of the missing species in main
#' clade and subclade or a vector c(missnumspec_M,missnumspec_S) with missing
#' species in main clade and subclade respectively.
#' @param methode The method used to solve the master equation, default is
#' 'analytical' which uses matrix exponentiation; alternatively numerical ODE
#' solvers can be used, such as 'odeint::runge_kutta_cash_karp54'. These were used in the
#' package before version 3.1.
#' @return The loglikelihood
#' @author Rampal S. Etienne & Bart Haegeman
#' @seealso \code{\link{dd_KI_ML}}, \code{\link{dd_loglik}}
#' \code{\link{dd_SR_loglik}}
#' @references - Etienne, R.S. et al. 2012, Proc. Roy. Soc. B 279: 1300-1309,
#' doi: 10.1098/rspb.2011.1439 \cr - Etienne, R.S. & B. Haegeman 2012. Am. Nat.
#' 180: E75-E89, doi: 10.1086/667574
#' @keywords models
#' @examples
#'
#' pars1 = c(0.25,0.12,25.51,1.0,0.16,8.61,9.8)
#' pars2 = c(200,1,0,18.8,1,2)
#' missnumspec = 0
#' brtsM = c(25.2,24.6,24.0,22.5,21.7,20.4,19.9,19.7,18.8,17.1,15.8,11.8,9.7,8.9,5.7,5.2)
#' brtsS = c(9.6,8.6,7.4,4.9,2.5)
#' dd_KI_loglik(pars1,pars2,brtsM,brtsS,missnumspec)
#'
#' @export dd_KI_loglik
dd_KI_loglik <- function(pars1,
pars2,
brtsM,
brtsS,
missnumspec,
methode = 'odeint::runge_kutta_cash_karp54')
{
if(length(pars2) == 4)
{
pars2[5] = 0
pars2[6] = 2
pars2[7] = 1
}
if(is.na(pars2[7]))
{
pars2[7] = 0
}
#pars2 <- as.data.frame(rbind(pars2))
names(pars2) = c('lx','ddep','cond','t_split','verbose','soc','corr')
tinn = -abs(pars1[7])
soc = pars2[6]
#pars2 <- as.vector(pars2)
# order branching times
brts = -sort(abs(c(brtsM,brtsS)),decreasing = TRUE)
if(sum(brts == 0) == 0)
{
brts[length(brts) + 1] = 0
}
S = length(brts) + (soc - 2)
brtsM = -sort(abs(brtsM),decreasing = TRUE)
brtsS = -sort(abs(brtsS),decreasing = TRUE)
tinn <- -abs(pars1[7])
if(!any(brtsM == 0))
{
brtsM = c(brtsM,0)
}
if(!any(brtsS == 0))
{
brtsS = c(brtsS,0)
}
if(!any(brtsS == tinn))
{
brtsS = c(tinn,brtsS)
}
# avoid coincidence of branching time and key innovation time
if(sum(abs(brtsM - tinn) < 1E-14) == 1) { tinn = tinn - 1E-8 }
ka = sum(brtsM < tinn)
brtsMtinn <- sort(c(tinn,brtsM),decreasing = FALSE)
kvec <- c(soc:(soc + ka - 1),
(soc + ka - 2):(soc + length(brtsMtinn) - 4),
soc + length(brtsMtinn) - 4)
brts_kM <- rbind(brtsMtinn,kvec)
kvec <- c(1:(length(brtsS)-1),length(brtsS)-1)
brts_kS <- rbind(brtsS,kvec)
brts_k_list <- list(brts_kM,brts_kS)
pars1_list <- list(pars1[1:3],pars1[4:6])
missnumspec_list <- create_missnumspec_list(missnumspec)
are_pars_impossible <- check_for_impossible_pars_KI(S_list = create_S_list(brts_k_list,soc),
missnumspec_list = missnumspec_list,
pars1_list = pars1_list,
pars2 = pars2)
if(are_pars_impossible || tinn >= 0 || tinn <= min(brtsM))
{
loglik <- -Inf
return(loglik)
}
loglik <- dd_multiple_KI_loglik(pars1_list = pars1_list,
pars2 = pars2,
brts_k_list = brts_k_list,
missnumspec_list = missnumspec_list,
reltol = 1e-14,
abstol = 1e-16,
methode = methode)
if(pars2[5])
{
s1 = sprintf('Parameters: %f %f %f %f %f %f %f, ',pars1[1],pars1[2],pars1[3],pars1[4],pars1[5],pars1[6],pars1[7])
s2 = sprintf('Loglikelihood: %f',loglik)
cat(s1,s2,"\n",sep = "")
utils::flush.console()
}
loglik = as.numeric(loglik)
if(is.nan(loglik) || is.na(loglik))
{
loglik = -Inf
}
return(loglik)
}
shift <- function(probs,pars2,k1)
{
lp <- length(probs)
if(pars2[7] == 1)
{
probs <- probs * 1/(k1 + (0:(lp - 1)))
} else
if(pars2[7] == 1.5)
{
probs <- probs * k1/(k1 + (0:(lp - 1)))
} else
if(pars2[7] == 2)
{
probs <- probs * (0:(lp - 1))/(k1 + (0:(lp - 1)))
} else
if(pars2[7] == 3)
{
probs[1:(lp - 1)] <- probs[2:lp] * (1:(lp - 1))/(k1 + (1:(lp - 1)))
probs[lp] <- 0
}
return(probs)
}
dd_KI_loglik_partial <- function(brts_k,
m,
cp,
pars1,
pars2,
lx,
reltol = 1e-14,
abstol = 1e-16,
methode,
log = TRUE)
{
ddep <- pars2[2]
verbose <- pars2[5]
for(i in 2:length(brts_k[1,]))
{
loglik <- cp[[1]]
probs <- cp[[2]]
t1 <- brts_k[1,i - 1]
t2 <- brts_k[1,i]
k1 <- brts_k[2,i - 1]
k2 <- brts_k[2,i]
y <- dd_integrate(probs,
c(t1,t2),
'dd_loglik_rhs',
c(pars1,k1,ddep),
rtol = reltol,
atol = abstol,
method = methode)
probs <- y[2,2:(lx + 1)]
if(k2 > k1)
{
probs <- flavec(ddep,
pars1[1],
pars1[2],
pars1[3],
r = 0,
lx,
k1) * probs # speciation event
} else if(k2 <= k1 && t2 != 0)
{
probs <- shift(probs = probs,
pars2 = pars2,
k1 = k1)
}
cp <- check_probs(loglik = loglik,
probs = probs,
verbose = verbose)
}
loglik <- cp[[1]]
probs <- cp[[2]]
if(log == TRUE)
{
loglik <- loglik + log(probs[1 + m])
} else
{
loglik <- exp(loglik) * probs[1 + m]
}
return(loglik)
}
convolve_logliks <- function(missnumspec_list,loglik_list)
{
if(length(missnumspec_list[[1]]) == 1)
{
loglik <- sum(unlist(loglik_list))
} else
{
loglik <- 0
llist <- length(loglik_list)
loglikdelta <- list()
convolve_lik <- 1
for(i in 1:llist)
{
loglikmax <- max(loglik_list[[i]])
loglikdelta[[i]] <- loglik_list[[i]] - loglikmax
loglik <- loglik + loglikmax
convolve_lik <- conv(convolve_lik,exp(loglikdelta[[i]]))
}
loglik <- loglik + log(convolve_lik[max(missnumspec_list[[1]]) + 1])
}
if(is.nan(loglik) || is.na(loglik))
{
loglik <- -Inf
}
return(loglik)
}
dd_KI_logliknorm <- function(brts_k_list,
pars1_list,
loglik,
cond,
ddep,
lx_list,
reltol = 1e-14,
abstol = 1e-16,
methode)
{
if(cond == 0 || loglik == -Inf)
{
logliknorm = 0
} else {
if(length(brts_k_list) > 2 && cond == 1)
{
stop('Conditioning on survival not implemented for more than 1 shift')
}
if(cond == 2 || cond == 3)
{
stop('Conditioning on number of species not implemented')
}
tcrown <- min(brts_k_list[[1]][1,])
ln <- length(brts_k_list[[1]][1,])
kinn <- 1 + min(which(brts_k_list[[1]][2,2:ln] - brts_k_list[[1]][2,1:(ln - 1)] < 0))
tinn <- brts_k_list[[1]][1,kinn]
tpres <- 0
# COMPUTE NORMALIZATION
# compute survival probability of clade S
lx = lx_list[[2]]
nx = -1:lx
lambdamu_nk <- lambdamu(nx,pars = c(pars1_list[[2]],0),ddep = ddep)
lavec <- lambdamu_nk[[1]]
muvec <- lambdamu_nk[[2]]
probs = rep(0,lx) # probs[1] = extinction probability
probs[2] = 1 # clade S starts with one species
if(methode != 'analytical')
{
m1 = lavec[1:lx] * nx[1:lx]
m2 = muvec[3:(lx + 2)] * nx[3:(lx + 2)]
m3 = (lavec[2:(lx + 1)] + muvec[2:(lx + 1)]) * nx[2:(lx + 1)]
if (startsWith(methode, 'odeint::')) {
probs = dd_logliknorm1_odeint(probs, c(tinn,tpres), c(m1,m2,m3), abstol, reltol, methode)
}
else {
y = deSolve::ode(probs,c(tinn,tpres),dd_logliknorm_rhs1,c(m1,m2,m3),rtol = reltol,atol = abstol,method = methode)
probs = y[2,2:(lx+1)]
}
} else
{
probs = dd_loglik_M(pars1_list[[2]],lx,0,ddep,tt = abs(tpres - tinn),probs)
}
PS = 1 - probs[1]
# compute survival probability of clade M
lx = lx_list[[1]]
n <- -1:lx
nx1 = rep(n,lx + 2)
dim(nx1) = c(lx + 2,lx + 2) # row index = number of species in first group
nx2 = t(nx1) # column index = number of species in second group
nxt = nx1 + nx2
lambdamu_n <- lambdamu2(nxt,pars1_list[[1]],ddep)
lavec <- lambdamu_n[[1]]
muvec <- lambdamu_n[[2]]
probs = matrix(0,lx,lx)
# probs[1,1] = probability of extinction of both lineages
# sum(probs[1:lx,1]) = probability of extinction of second lineage
probs[2,2] = 1 # clade M starts with two species
# STEP 1: integrate from tcrown to tinn
dim(probs) = c(lx,lx)
if(methode != 'analytical')
{
m1 = lavec[1:lx,2:(lx+1)] * nx1[1:lx,2:(lx+1)]
m2 = muvec[3:(lx+2),2:(lx+1)] * nx1[3:(lx+2),2:(lx+1)]
ma = lavec[2:(lx+1),2:(lx+1)] + muvec[2:(lx+1),2:(lx+1)]
m3 = ma * nx1[2:(lx+1),2:(lx+1)]
m4 = lavec[2:(lx+1),1:lx] * nx2[2:(lx+1),1:lx]
m5 = muvec[2:(lx+1),3:(lx+2)] * nx2[2:(lx+1),3:(lx+2)]
m6 = ma * nx2[2:(lx+1),2:(lx+1)]
if (startsWith(methode, "odeint::")) {
probs = dd_logliknorm2_odeint(probs, c(tcrown,tinn), list(m1,m2,m3,m4,m5,m6), reltol, abstol, methode)
}
else {
dim(probs) = c(lx*lx,1)
y = deSolve::ode(probs,c(tcrown,tinn),dd_logliknorm_rhs2,list(m1,m2,m3,m4,m5,m6),rtol = reltol,atol = abstol, method = "ode45")
probs = y[2,2:(lx * lx + 1)]
}
} else
{
probs = dd_loglik_M2(pars = pars1_list[[1]],lx = lx,ddep = ddep,tt = abs(tinn - tcrown),p = probs)
}
dim(probs) = c(lx,lx)
probs[1,1:lx] = 0
probs[1:lx,1] = 0
# STEP 2: transformation at tinn
nx1a = nx1[2:(lx+1),2:(lx+1)]
nx2a = nx2[2:(lx+1),2:(lx+1)]
probs = probs * nx1a/(nx1a+nx2a)
probs = rbind(probs[2:lx,1:lx], rep(0,lx))
# STEP 3: integrate from tinn to tpres
if(methode != 'analytical')
{
if (startsWith(methode, "odeint::")) {
probs = dd_logliknorm2_odeint(probs, c(tinn, tpres), list(m1,m2,m3,m4,m5,m6), reltol, abstol, 'odeint::runge_kutta_fehlberg78')
}
else {
dim(probs) = c(lx*lx,1)
y = deSolve::ode(probs,c(tinn, tpres),dd_logliknorm_rhs2,list(m1,m2,m3,m4,m5,m6),rtol = reltol,atol = abstol, method = "ode45")
probs = y[2,2:(lx * lx + 1)]
dim(probs) = c(lx,lx)
}
} else
{
probs = dd_loglik_M2(pars = pars1_list[[1]],lx = lx,ddep = ddep,tt = abs(tpres - tinn),p = probs)
}
dim(probs) = c(lx,lx)
PM12 = sum(probs[2:lx,2:lx])
PM2 = sum(probs[1,2:lx])
#print(log(2) + log(PM12 + PS * PM2))
#print(log(2) + (log(PM12 + PM2) + log(PS)))
#print(log(2) + (log(PM12) + log(PS)))
logliknorm = log(2) + (cond == 1) * log(PM12 + PS * PM2) +
(cond == 4) * (log(PM12 + PM2) + log(PS)) +
(cond == 5) * (log(PM12) + log(PS))
}
return(logliknorm)
}
create_lx_list <- function(lmax,
ddep,
missnumspec_list,
pars1_list)
{
lx_list <- list()
for(i in 1:length(pars1_list))
{
if(ddep == 1)
{
lx_list[[i]] <- min(max(1 + max(missnumspec_list[[i]]),1 + ceiling(pars1_list[[i]][1]/(pars1_list[[i]][1] - pars1_list[[i]][2]) * pars1_list[[i]][3])),ceiling(lmax))
} else {
if(ddep == 1.3)
{
lx_list[[i]] <- min(max(1 + max(missnumspec_list[[i]]),1 + ceiling(pars1_list[[i]][3])),ceiling(lmax))
} else {
lx_list[[i]] <- round(lmax)
}
}
}
return(lx_list)
}
create_missnumspec_list <- function(m)
{
if(length(m) == 1)
{
missnumspec_list <- list(0:m,0:m)
} else
{
missnumspec_list <- as.list(m)
}
return(missnumspec_list)
}
create_S_list <- function(brts_k_list,soc)
{
S_list <- list()
for(i in 1:length(brts_k_list))
{
S_list[[i]] <- brts_k_list[[i]][2,length(brts_k_list[[i]][2,])]
names(S_list[[i]]) <- NULL
}
return(S_list)
}
check_for_impossible_pars_KI <- function(S_list,
missnumspec_list,
pars1_list,
pars2)
{
ddep <- pars2[2]
loglik <- 0
sumK <- 0
sumSm <- max(missnumspec_list[[1]])
for(i in 1:length(S_list))
{
are_pars_impossible <- check_for_impossible_pars(pars1_list[[i]],ddep)
if(are_pars_impossible)
{
return(are_pars_impossible)
}
m <- 0
if(length(missnumspec_list[[i]]) == 1)
{
m <- missnumspec_list[[i]]
}
are_pars_impossible <-
((ddep == 1 && ceiling(pars1_list[[i]][1]/(pars1_list[[i]][1] - pars1_list[[i]][2]) * pars1_list[[i]][3]) < S_list[[i]] + m) ||
(ddep == 1.3 && ceiling(pars1_list[[i]][3]) < S_list[[i]] + m))
if(are_pars_impossible)
{
return(are_pars_impossible)
}
sumK <- sumK + ceiling(pars1_list[[i]][3])
sumSm <- sumSm + S_list[[i]]
}
are_pars_impossible <- (ddep == 1.3 && sumK < sumSm)
return(are_pars_impossible)
}
check_for_impossible_pars <- function(pars1,
ddep)
{
are_pars_impossible <- min(pars1) < 0
if(are_pars_impossible)
{
cat('One or more parameters are negative.')
} else
{
are_pars_impossible <-
(pars1[2] == 0 && (ddep == 2 || ddep == 2.1 || ddep == 2.2)) ||
(pars1[1] == 0 && (ddep == 4 || ddep == 4.1)) ||
(pars1[1] <= pars1[2] && (ddep == 1 || ddep == 2 || ddep == 2.1 || ddep == 2.2 || ddep == 3 || ddep == 4 || ddep == 4.1))
if(are_pars_impossible)
{
cat('These parameters are incompatible with the chosen model')
}
}
return(are_pars_impossible)
}
#' Loglikelihood for diversity-dependent diversification models with multiple
#' decoupling (rate shift) events
#'
#' This function computes loglikelihood of a diversity-dependent
#' diversification model for a given set of branching times and parameter
#' values where the diversity-dependent dynamics of subclades decouple from
#' the dynamics of main clades, potentially accompanied by a
#' shift in parameters.
#' @inheritParams dd_KI_loglik
#' @param pars1_list list of paramater sets one for each rate regime (subclade).
#' The parameters are: lambda (speciation rate), mu (extinction rate), and K
#' (clade-level carrying capacity).
#' @param brts_k_list list of matrices, one for each rate regime (subclade). Each
#' matrix has in the first row the branching times including the shift/decoupling time
#' and the present time (0) in negative time (i.e. 10 mya = -10). In the second row
#' it has the number of lineages, i.e. starting at 2 for a phylogeny with a crown
#' and increasing by one at each branching time and decreasing by one at each
#' decoupling/shift time. The last element is the same as the second last.
#' @param missnumspec_list list containing the number of missing species for each clade.
#' If only a single number m of missing species is known for the entire phylogeny, then each
#' element of the list should be 0:m. One can also create this from m using the function
#' create_missnumspec_list
#' @param reltol relative tolerance in integration of the ODE system, default at 1e-14
#' @param abstol tolerance tolerance in integration of the ODE system, default at 1e-16
dd_multiple_KI_loglik <- function(pars1_list,
pars2,
brts_k_list,
missnumspec_list,
reltol = 1e-14,
abstol = 1e-16,
methode = 'odeint::runge_kutta_cash_karp54')
{
lx_list <- create_lx_list(lmax = pars2[1],
ddep = pars2[2],
missnumspec_list = missnumspec_list,
pars1_list = pars1_list)
llist <- length(pars1_list)
if(llist != length(brts_k_list) || llist != length(missnumspec_list))
{
stop('The length of the lists of parameters, branching times and missing species is not the same.')
}
if(llist < 2)
{
stop('The length of the lists of parameters must be 2 or larger.')
}
loglik_list <- as.list(rep(0,llist))
for(i in 1:llist)
{
probs <- rep(0,lx_list[[i]])
probs[1] <- 1
loglik_list[[i]] <- dd_KI_loglik_partial(brts_k = brts_k_list[[i]],
m = missnumspec_list[[i]],
cp = list(0,probs),
pars1 = pars1_list[[i]],
pars2 = pars2,
lx = lx_list[[i]],
reltol = reltol,
abstol = abstol,
methode = methode)
}
loglik <- convolve_logliks(missnumspec_list,loglik_list)
if(pars2[3] != 5)
{
logliknorm <- dd_KI_logliknorm(brts_k_list = brts_k_list,
pars1_list = pars1_list,
loglik = loglik,
cond = pars2[3],
ddep = pars2[2],
lx_list = lx_list,
reltol = reltol,
abstol = abstol,
methode = methode)
} else
{
pars2[7] <- 3
logliknorm <- dd_multiple_KI_logliknorm(brts_k_list = brts_k_list,
pars1_list = pars1_list,
pars2 = pars2,
loglik = loglik,
lx_list = lx_list,
reltol = 1e-14,
abstol = 1e-16,
methode = methode)
}
loglik <- loglik - logliknorm
S_list <- create_S_list(brts_k_list = brts_k_list,
soc = pars2[6])
Sv <- unlist(S_list)
m <- unlist(missnumspec_list)
if(length(missnumspec_list[[1]]) > 1)
{
m <- max(missnumspec_list[[1]])
Sv <- sum(Sv)
}
loglik <- loglik - sum(lgamma(Sv + m + 1) - lgamma(Sv + 1) - lgamma(m + 1))
return(loglik)
}
dd_multiple_KI_logliknorm <- function(brts_k_list,
pars1_list,
pars2,
loglik,
lx_list,
reltol = 1e-14,
abstol = 1e-16,
methode)
{
if(pars2[3] == 0 || loglik == -Inf)
{
logliknorm <- 0
} else
{
llist <- length(pars1_list)
loglik_list <- as.list(rep(0,llist))
logliknorm <- 0
for(i in 1:llist)
{
probs <- rep(0,lx_list[[i]])
probs[1] <- 1
tstart <- min(brts_k_list[[i]][1,])
ln <- length(brts_k_list[[i]][1,])
kinn <- 1 + (which(brts_k_list[[i]][2,2:ln] - brts_k_list[[i]][2,1:(ln - 1)] < 0))
tinn <- brts_k_list[[i]][1,kinn]
tpres <- 0
brts_k <- rbind(c(tstart,tinn,tpres),rep(brts_k_list[[i]][2,1],length(tinn) + 2))
lp <- length(probs)
if(brts_k_list[[i]][2,1] == 1)
{
aux <- 1:lp
}
if(brts_k_list[[i]][2,1] == 2)
{
aux <- (2:(lp + 1)) * (3:(lp + 2))/6
}
loglik_list[[i]] <- dd_KI_loglik_partial(brts_k = brts_k,
m = 0:(lp - 1),
cp = list(0,probs),
pars1 = pars1_list[[i]],
pars2 = pars2,
lx = lx_list[[i]],
reltol = reltol,
abstol = abstol,
methode = methode,
log = FALSE)/aux
logliknorm <- logliknorm + log(sum(loglik_list[[i]]))
}
}
return(logliknorm)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.