# Stochastic logistic model
# with variation in R0 & K possible
N0=50; T=15
muR0=0.2; sdR0=0
muK=100; sdK=0
sdObs=0;
N=matrix(NA,T+1,1)
N[1]=N0
Nobs=matrix(NA,T+1,1)
Nobs[1]=N0+rnorm(1,0,sdObs)

for (t in 1:T){
	R0=rnorm(1,muR0,sdR0)
	K =rnorm(1,muK,sdK) 
	N[t+1]=N[t]*(1+R0*(1-N[t]/K))
	Nobs[t+1]=N[t+1]+rnorm(1,0,sdObs)
  }

Nmax=max(max(Nobs),max(N))
Nmin=min(min(Nobs),min(N))

plot(1:(T+1),N,"o",ylim=c(0,130),xlab='Year',
	main=c('avg N =', round(mean(N[10:T]),2)))
abline(h=muK,col='red',lwd=3)
points(1:(T+1),Nobs,"o",col='blue')