#    Calculates by simulation the probability that a population
#    undergoing density-dependent, theta-logistic growth will fall 
#    below the quasi-extinction threshold Nx at or before time tmax
#    Assumes that a population starts at its carrying capacity (K)

############# Simulation input parameters ##########################

r <- 0.5967		# intrinsic rate of increase
theta <- 1		# nonlinearity in density dependence
sigma2 <- 0.0745	# environmental variance of r
maxN <- 30000  		# maximum initial population size (Nc) to simulate
interval <- 100		# sample possible Nc values how often
probcast <- 0.2		# per year probability of catastrophe
alpha <- 0.5		# the shape1 (alpha) parameter from a Beta distribution characterizing the distribution of fish-kill severities
beta <- 2.8		# the shape2 (beta) parameter from a Beta distribution characterizing the distribution of fish-kill severities
Nx <- 50		# quasi-extinction threshold
tmax <- 100		# Number of years being simulated
NumReps <- 10000	# number of replicate population trajectories per Nc value


############# Population Growth Simulation #########################

sigma <- sqrt(sigma2)
inter <- interval*(c(1:(maxN/interval)))
M <- matrix(NA,nrow=length(inter),ncol=3)
colnames(M) <- c("K","Prob.persist","Prob.extinct")
for (j in 1:length(inter)) {
N <- matrix(NA,nrow=tmax,ncol=NumReps)
K <- j*interval
Nc <- K
N[1,] <- Nc
for (i in 1:NumReps) {								
  for (t in 2:tmax) {									
    catastrophe <- 1 - (rbinom(1*1, size=1, prob=probcast)*rbeta(1*1, shape1=alpha, shape2=beta))
    N[t,i] <- N[(t-1),i]*exp(r*(1-((N[(t-1),i]/K)^theta))+(sigma*rnorm(1,mean=0,sd=1)))*catastrophe					
  }
}
M[j,1] <- K
M[j,2] <- 1-((NumReps-length(which(N[tmax,]>0)))/NumReps)
M[j,3] <- 1-M[j,2]
}

############# Estimate of MVP ######################################

MVP <- interval*min(which(M[,2]>0.94999))
MVP
write.table(M,"MVP.txt")

############# Plot of extinction probability #######################

plot(M[,1],M[,3], 'l', col = 'blue', ylim=c(0,1), xlab='Initial population size', ylab='Probability of extinction')