In franciscorichter/dmea02: Diversification and Macro-Evolution Analysis

knitr::opts_chunk$set(echo = TRUE)
library(subplex)

S1. Simulate lighting data

sim.light <- function(T,lambda=0.8,p=0.9){
  wt = NULL
  tau = NULL
  obs = NULL
  ot=0
  i = 1
  j=1
  while(sum(wt)<T){
   wt[i] = rexp(1,rate = lambda)
   tau[i] = rbinom(1,1,prob=p)
   ot = ot+wt[i]
   if(tau[i]==1){
     obs[j] = ot
     ot = 0
     j=j+1
   }
   i=i+1
  }
  d = length(wt)
  wt = wt[-d]
  tau = tau[-d]
  #wt[d] = T-sum(wt)
  return(list(wt=wt,tau=tau,obs=obs))
}


l1 = sim.light(50,p=0.5)
l1

S.1.1 Calculate MLE

llik.light <- function(light,par){
  lambda = par[1]
  p = par[2]
  exp = dexp(x=light$wt,rate=lambda,log=T)
  bin = dbinom(x=sum(light$tau),size=length(light$tau),prob=p,log=TRUE)
  llik = sum(exp)+bin
  if(p>1 | p<0){
    llik = -Inf
  }
  return(llik)
}
mle.light <- function(light,init_par = c(2,1)){
  p = sum(light$tau)/length(light$tau)
  lambda = 1/mean(light$wt)
  #optim(par=c(0.5,0.5),fn=llik.light,lower=c(0,0),light=light, upper = c(5,1),method="L-BFGS-B")
  return(list(lambda=lambda,p=p))
}
l1
mle.light(l1,init_par=c(1,0.5))

n=100
L = vector(mode='numeric',length = n)
P = vector(mode='numeric',length = n)
for(i in 1:n){
  li = sim.light(T=50)
  ml = mle.light(light = li)
  L[i] = ml$lambda
  P[i] = ml$p
}

summary(L)
summary(P)

Incomplete data

Option 1. The current way

rec.light <- function(obslight,lambda,p){
  dt = obslight
  for(i in 1:length(dt)){
    t = dt[i]

  }
}

Option 2. topoloy first

sim.miss.light <- function(obslight,lambda,p){
  dt = obslight
  j = 1
  wt = NULL
  top = NULL
  error = NULL
  timesuggestion = NULL
  for(i in 1:length(dt)){

    to = m = 0
    while(to == 0){
      to = rbinom(1,1,p)
      top[j] = to
      j = j+1
      if(to == 0) m <- m+1
    }
    exps = rexp(n=m+1,rate=lambda)
    swt = sum(exps)
    error[i] = abs(obslight[i]-swt)/(m+1)
    timesuggestion[i] = obslight[i]<swt
    wt = c(wt,exps)
  }
  light = list(wt=wt,tau=top,obs=obslight,error=error,errorsum=summary(error),timesuggestion=timesuggestion)
  return(light)
}


s = sim.light(50,p=0.5)
s
rec = sim.miss.light(obslight = s$obs,lambda=0.8,p=0.5)
rec
n=10000
min = vector(mode='numeric',length=n)
max = vector(mode='numeric',length=n) 
med = vector(mode='numeric',length=n)  
lambda = vector(mode='numeric',length=n)  
p = vector(mode='numeric',length=n)  
for(i in 1:n){
  rec = sim.miss.light(obslight = s$obs,lambda=0.8,p=0.5)
  min[i] = rec$errorsum[1]
  med[i] = rec$errorsum[3]
  max[i] = rec$errorsum[6]
  mle = mle.light(light = rec)
  lambda[i] = mle$lambda
  p[i] = mle$p
}
library(ggplot2)
qplot(med,geom='histogram',binwidth=0.1)
summary(lambda)
summary(p)


###


n=10000
min = vector(mode='numeric',length=n)
max = vector(mode='numeric',length=n) 
med = vector(mode='numeric',length=n)  
lambda = vector(mode='numeric',length=n)  
p = vector(mode='numeric',length=n)  
for(i in 1:n){
  rec = sim.miss.light(obslight = s$obs,lambda=0.4,p=0.5)
  min[i] = rec$errorsum[1]
  med[i] = rec$errorsum[3]
  max[i] = rec$errorsum[6]
  mle = mle.light(light = rec)
  lambda[i] = mle$lambda
  p[i] = mle$p
}
library(ggplot2)
qplot(med,geom='histogram',binwidth=0.1)
summary(lambda)
summary(p)



###


n=10000
min = vector(mode='numeric',length=n)
max = vector(mode='numeric',length=n) 
med = vector(mode='numeric',length=n)  
lambda = vector(mode='numeric',length=n)  
p = vector(mode='numeric',length=n)  
for(i in 1:n){
  rec = sim.miss.light(obslight = s$obs,lambda=0.8,p=0.2)
  min[i] = rec$errorsum[1]
  med[i] = rec$errorsum[3]
  max[i] = rec$errorsum[6]
  mle = mle.light(light = rec)
  lambda[i] = mle$lambda
  p[i] = mle$p
}
library(ggplot2)
qplot(med,geom='histogram',binwidth=0.1)
summary(lambda)
summary(p)

Option 3. Times first

sim.miss.light <- function(obslight,lambda,p){
  dt = obslight
  j = 1
  wt = NULL
  top = NULL
  error = NULL
  for(i in 1:length(dt)){

    to = m = 0
    while(sum(exp)<dt[i]){
      to = rexp
      top[j] = to
      j = j+1
      if(to == 0) m <- m+1
    }
    exps = rexp(n=m+1,rate=lambda)
    swt = sum(exps)
    error[i] = abs(obslight[i]-swt)/(m+1)
    timesuggestion[i] = obslight[i]<swt
    wt = c(wt,exps)
  }
  light = list(wt=wt,tau=top,obs=obslight,error=error,errorsum=summary(error),timesuggestion=timesuggestion)
  return(light)
}


s = sim.light(50,p=0.5)
s
rec = sim.miss.light(obslight = s$obs,lambda=0.8,p=0.5)
rec