model <- function(t, state, parms) { with(as.list(c(state,parms)), { R0 <- b/d k1 <- K/(2*(1-1/R0)) k2 <- K/(R0-1) k3 <- K/sqrt(R0-1) k4 <- K*log(2)/log(R0) f1 <- 1 - N1/(2*k1) f2 <- 1/(1+N2/k2) f3 <- 1/(1+(N3/k3)^2) f4 <- exp(-log(2)*N4/k4) dN1 <- b*N1*f1 - d*N1 dN2 <- b*N2*f2 - d*N2 dN3 <- b*N3*f3 - d*N3 dN4 <- b*N4*f4 - d*N4 return(list(c(dN1,dN2,dN3,dN4))) }) } p <- c(b=0.2,d=0.1,K=1) s <- c(N1=0.01,N2=0.01,N3=0.01,N4=0.01) run() #How well do these 4 functions buffer noise? s <- c(N1=1,N2=1,N3=1,N4=1) after="state<-state+rnorm(1,0,0.1);state<-ifelse(state<0,0,state)" data <- run(500,after=after,table=T) apply(data,2,sd) after="parms[3]<-abs(rnorm(1,1,0.5))" data <- run(500,after=after,table=T) apply(data,2,sd)