test omega for mgdrive

getOmega.R

get_omega1 <- function(mu,lifespanReduction){
  lifespan <- (1/mu) * lifespanReduction
  get_omega_f <- function(omega,mu,lifespan){
    (1/(1-omega+(omega*mu)))-lifespan
  }
  return(uniroot(f=get_omega_f,interval=c(0,1),lifespan=lifespan,mu=mu,maxiter=1e4)$root)
}

# 10-25-2018
get_omega <- function(mu,lifespan){
  # function to find root
  fn <- function(mu,lifespan,omega){
    (1 / (1 - ((1 - mu) * omega))) - lifespan
  }
  bound <- 1 / (1 - mu)
  return(uniroot(f=fn,interval=c(0,bound),lifespan=lifespan,mu=mu,maxiter=1e4)$root)
}

popSim <- function(init,t,mu,omega){
  pop <- rep(0,t)
  pop[1] <- init
  for(i in 2:t){
    pop[i] <- pop[i-1]*((1-mu)*omega)
  }
  return(pop)
}

cdfG <- function(p,x){
  1-((1-p)^x)
}

mu <- 1/10

tmax <- 50
out <- popSim(init = 100,t = tmax,mu = mu,omega = 1)
plot(x = 1:tmax,y = out,type = "l")
grid()

out[25]/out[1]
1-cdfG(mu,24)

omega <- get_omega1(mu = mu,lifespanReduction = 0.75)
# omega <- MGDrivE::get_omega(mu = mu,lifespanReduction = 0.75)
out1 <- popSim(init = 100,t = tmax,mu = mu,omega = omega)

lines(x = 1:tmax,y = out1,col="red")

out1[25]/out1[1]
1-cdfG((1-omega+(omega*mu)),24)