############################### # linear stochastic simulator # ############################### # a simulator for the model # N(t+1)=r(t)*N(t) # where r(t) is given by an empirical distribution # INPUT: # r the empirical distribution of r values (choosen with equal likelihood) # reps the number of simulations to run # T the length of each simulation # N0 initial density # plottype 0 means none, 1 means population trajectories on regular scale, 2 means log densities and histogram linear=function(r=c(2,0.2),reps=1,T=50,N0=100,plottype=1){ N=matrix(N0,T,reps) for(i in 1:(T-1)){ N[i+1,]=sample(r,reps,replace=TRUE)*N[i,]} if(plottype==1){par(cex.lab=1.5,cex.axis=1.5); matplot(N,type="l",pch=21,bg="red",lwd=3,xlab="time",ylab="N")} if(plottype==2){par(cex.lab=1.5,cex.axis=1.5,mfrow=c(1,2),mar=c(4,4,1,1)); matplot(log(N),type="l",xlab="time",ylab="log(density)") abline(a=log(N0),b=mean(log(r)),lwd=4,lty=2) hist(log(N[T,]),40,col="red",yaxt="n",ylab="",xlab="log(density)",main="") } } ######################## ## Bethedging Function # ######################## # this function determines the germination fraction g that maximizes # the expected log fitness log(gY+(1-g)s) where # Y is a discrete random variable taking on the values Y=c(....) # with probabilities p=c(....) and s is probability of a seed suviving # INPUT # Y list of yield values # p list of probabilities of yield values # s seed survival # plotit - TRUE creates a plot, FALSE supresses the plot # OUTPUT # say out=bethedging(....) then # out\$maximum is the g value that maximizes the log-fitness # out\$objective is the maximal log-fitness hedge=function(Y=c(100,0),p=c(0.5,0.5),s=0.9,plotit=TRUE){ f=function(g)sum(p*log(g*Y+(1-g)*s)); temp=optimize(f,interval=c(0,1),maximum=TRUE); if(plotit){ g=seq(0,0.999999,length=100) out=sapply(g,f) par(cex.lab=1.25,cex.axis=1.25,mfrow=c(1,1),mar=c(4.5,4.5,1,1)) plot(g,out,type="l",lwd=3,col="red",xlab="g",ylab="log(fitness)") abline(h=0,lty=2) abline(v=temp\$maximum,lty=2,col="blue") abline(h=temp\$objective,lty=2,col="blue")} return(temp);} ################################### ## Stochastic Beverton-Holt model # ################################### # a simulator for the model # N(t+1)=r(t)*N(t)/(1+a*N[t]) # where r(t) is given by an empirical distribution # and a is the strength of density dependence # INPUT: # r the empirical distribution of r values (choosen with equal likelihood) # a strength of intraspecific competition # reps the number of simulations to run # T the length of each simulation # N0 initial density # plottype 0 means none # 1 means population trajectories on regular scale # 2 means log densities and histogram of all trajectories at time T # 3 means densties and historgram of first trajectrom from T/2 to T nonlinear=function(r=c(1.2,0.95),reps=1,T=200,N0=1,plottype=1){ N=matrix(N0,T,reps) for(i in 1:(T-1)){ N[i+1,]=sample(r,reps,replace=TRUE)*N[i,]/(1+0.001*N[i,])} if(plottype==1){par(cex.lab=1.5,cex.axis=1.5); matplot(N,type="l",pch=21,bg="red",lwd=3,xlab="time",ylab="N")} if(plottype==2){par(cex.lab=1.5,cex.axis=1.5,mfrow=c(1,2),mar=c(4,4,1,1)); matplot(log(N),type="l",xlab="time",ylab="log(density)") hist(log(N[T,]),40,col="red",yaxt="n",ylab="",xlab="log(density)",main="") } if(plottype==3){par(cex.lab=1.5,cex.axis=1.5,mfrow=c(1,2),mar=c(4,4,1,1)); matplot(N[,1],type="l",xlab="time",ylab="density") hist(N[(T/2):T,1],40,xlim=c(0,max(N[,1])),col="red",yaxt="n",ylab="",xlab="density",main="") } }