R/Simulation_functions.R

In franciscorichter/emphasisR: Expectation Maximization on PHylogenetic Applications with Simulations and Importance Sampling

sim_brts_bootstrap <- function(pars,model="rpd5c",ct,bootstrap_n=1){
  
  sim_tree = get(paste0("sim.tree_",model))
  SIMS = vector(mode = "list",length = bootstrap_n)
  sim = NULL
  for(i in 1:bootstrap_n){
    sim$brts = 0
    number_of_empty_trees = -1
    while(length(sim$brts) == 1){
      sim = sim_tree(pars = pars,ct=ct,soc=2)
      number_of_empty_trees = number_of_empty_trees + 1 
    }
    SIMS[[i]] = c(sim,list(number_of_empty_trees=number_of_empty_trees))

  }
    
  return(SIMS)
}

### simulation of trees 

sim.tree_rpd1 <- function(pars,ct,timeLimit=10){
  setTimeLimit(timeLimit)
  tree =  data.frame(brts=c(0,0),
                     to=c(1,1),
                     t_ext=c(Inf,Inf),
                     lambda=c(0,0), 
                     clade = c("a","b"))
  cbt = 0
  N1 = 1 
  N2 = 1
  mu = max(0,pars[1])
  while((cbt < ct)  &  (N1 >= 1) &  (N2 >= 1)){
    N = N1 + N2 
    lambda_ct = max(0,pars[2] + pars[3]*N)  # diversity dependence only 
    rate_max = (lambda_ct+mu)*N
    u1 = runif(1)
    next_event_time = cbt-log(x = u1)/rate_max
    
    if(next_event_time < ct){

        to = sample(c(1,0),size=1,prob=c(lambda_ct,mu)/(lambda_ct+mu))
        clade = sample(c("a","b"),size=1,prob=c(N1,N2)/(N1+N2))
        if(to == 1){
          if(clade=="a"){
            N1 = N1+1
          }else{
            N2 = N2+1
          }
        }else{
          if(clade=="a"){
            N1 = N1-1
          }else{
            N2 = N2-1
          } 
          if((N1 >= 1) &  (N2 >= 1)){
            ext_spec = sample(which(tree$to == 1 & tree$t_ext == Inf & tree$clade == clade),1)
            tree$t_ext[ext_spec] = next_event_time
          }
        }
        tree = rbind(tree,data.frame(brts=next_event_time,to=to,t_ext=Inf,lambda=lambda_ct, clade = clade))
    }
    cbt = next_event_time
  }

  if((N1 == 0) |  (N2 == 0)){
    tree = 0
    brts = 0
  }else{
    tree = rbind(tree,data.frame(brts=ct,to=1,t_ext=Inf,lambda=lambda_ct, clade=0))
    brts = tree$brts[is.infinite(tree$t_ext) & tree$to==1]
  }
  return(list(tree=tree,brts=brts))
}


sim.tree_rpd5 <- function(pars,ct,timeLimit=2){
  setTimeLimit(timeLimit)
  tree =  data.frame(brts=c(0,0),
                     to=c(1,1),
                     t_ext=c(Inf,Inf),
                     clade = c("a","b"))
  cbt = 0
  N1 = 1 
  N2 = 1
  mu = max(0,pars[1])
  P=0
  while((cbt < ct)  &  (N1 >= 1) &  (N2 >= 1)){
    next_bt = min(c(tree$brts[tree$brts>cbt],ct))
    N = N1+N2
    lambda_mx = max(0,pars[2] + pars[3]*N  +  ((P + N*(next_bt-cbt)-cbt)/N )*pars[4])
    rate_max = (lambda_mx+mu)*N 
    u1 = runif(1)
    next_event_time = cbt-log(x = u1)/rate_max
    P = P + N*(next_event_time-cbt)
    
    if(next_event_time < next_bt){
      u2 = runif(1)
      lambda_ct = max(0,pars[2] + pars[3]*N  +  ((P+N*(next_event_time-cbt) - next_event_time)/N)*pars[4])
      pt =( (lambda_ct+mu)*N )/rate_max  
      
      if(u2<pt){
        to = sample(c(1,0),size=1,prob=c(lambda_ct,mu)/(lambda_ct+mu))
        clade = sample(c("a","b"),size=1,prob=c(N1,N2)/(N1+N2))
        if(to == 1){
          if(clade=="a"){
            N1 = N1+1
          }else{
            N2 = N2+1
          }
        }else{
          if(clade=="a"){
            N1 = N1-1
          }else{
            N2 = N2-1
          }
          if((N1 >= 1) &  (N2 >= 1)){
            ext_spec = sample(which(tree$to == 1 & tree$t_ext == Inf & tree$clade == clade),1)
            tree$t_ext[ext_spec] = next_event_time
            removed_branch_length = tree$t_ext[ext_spec] - tree$brts[ext_spec]
            P = P - removed_branch_length
          }
        }
            
        tree = rbind(tree,data.frame(brts=next_event_time,to=to,t_ext=Inf,clade=clade))
        
      }
    }
    
    cbt = next_event_time
  }
  if((N1 == 0) |  (N2 == 0)){
    tree = 0
    brts = 0
  }else{
    tree = rbind(tree,data.frame(brts=ct,to=1,t_ext=Inf, clade=0))
    brts = tree$brts[is.infinite(tree$t_ext) & tree$to==1]
  }
  return(list(tree=tree,brts=brts))
}

sim.tree_rpd5c <- function(pars,ct,soc=2){
  tree = NULL
  cbt = 0 
  N = soc
  mu = max(0,pars[1])
  P = 0
  while((cbt < ct)  &  (N >= soc)){
    lambda_mx = max(0,pars[2] + pars[3]*N  +  ((P + N*(ct-cbt)-cbt)/N )*pars[4])
    rate_max = (lambda_mx+mu)*N 
    u1 = runif(1)
    next_event_time = cbt-log(x = u1)/rate_max
    P = P + N*(next_event_time-cbt)
    if(next_event_time < ct){
      u2 = runif(1)
      lambda_ct = max(0,pars[2] + pars[3]*N  +  ((P+N*(next_event_time-cbt) - next_event_time)/N)*pars[4])
      pt =( (lambda_ct+mu)*N )/rate_max  
      
      if(u2<pt){
        to = sample(c(1,0),size=1,prob=c(lambda_ct,mu)/(lambda_ct+mu))
        if(to == 1){
          N = N + 1
        }else{
          N = N - 1 
          if(N>=soc){
            ext_spec = sample(which(tree$to == 1 & tree$t_ext == Inf),1)
            tree$t_ext[ext_spec] = next_event_time
            removed_branch_length = tree$t_ext[ext_spec] - tree$brts[ext_spec]
            P = P - removed_branch_length
          }
          
        }
        tree = rbind(tree,data.frame(brts=next_event_time,to=to,t_ext=Inf,lambda=lambda_ct))
        
      }
    }
    cbt = next_event_time
  }
  if(N==(soc-1)){
    tree = 0
    brts = 0
  }else{
    lambda_ct = max(0,pars[2] + pars[3]*N  +  ((P+N*(ct-cbt) - ct)/N)*pars[4])
    tree = rbind(tree,data.frame(brts=ct,to=1,t_ext=Inf,lambda=lambda_ct))
    brts = tree$brts[is.infinite(tree$t_ext) & tree$to==1]
  }
  return(list(tree=tree,brts=brts))
}



#### simulation conditional to n extant species
 
sim_cond_tree_cr <- function(pars,ct){
  
  lambda = pars[2]
  mu = pars[1]
  
  reject = 0
  while(reject==0){
    brts=n=NULL
    N = 1
    cbt = 0 
    while(cbt<ct & N>0){
      sigma = N*(lambda+mu)
      cbt = cbt + rexp(1,sigma)
      if(cbt < ct){
        a = sample(c(-1,1),size=1,prob=c(mu,lambda)/(mu+lambda))      
        N = N+a
        n = c(n,N)
        brts = c(brts,cbt)
      }else{
        if(N==1){
          reject=1
        }
      }
    }
  }
  return(data.frame(brts=c(brts,ct),n=c(1,n)))
}



sim_cond_tree <- function(pars,ct){
  reject = 0
  while(reject==0){
    brts=n=NULL
    N = 1
    cbt = 0 
    while(cbt<ct & N>0){
      lambda = max(0,pars[2]+pars[3]*N)
      mu = max(0,pars[1])
      
      sigma = N*(lambda+mu)
      cbt = cbt + rexp(1,sigma)
      if(cbt < ct){
        a = sample(c(-1,1),size=1,prob=c(mu,lambda)/(mu+lambda))      
        N = N+a
        n = c(n,N)
        brts = c(brts,cbt)
      }else{
        if(N==1){
          reject=1
        }
      }
    }
  }
  return(data.frame(brts=c(brts,ct),n=c(1,n)))
}