Mer <- 0.016242; Mir <- 0.008238; maxdate <- 365*2; PI <- rep(0,maxdate) a <- 5.0; delta <- 14.0; eta <- 14.0; r <- 5.0; omega <- 0.0001 mu <- 0.2; tau <- 10.0; phi <- 5.0; lambda <- 14.0 simexec <- function(px,beta) { ZI <- rep(0,maxdate) Ns <- 400; Ne <- 200; Ni <- 200; Nr <- 200 N <- Ns+Ne+Ni+Nr M <- floor(N*(1-px*beta)*r*a*((1-mu)^phi)/mu) Me <- floor(M*Mer); Mi <- floor(M*Mir); Ms <- M - Me - Mi for (j in 1:(maxdate-1)) { ZI[j] <- Ni/N IB <- rbinom(1,rbinom(1,Ns,(1-px*beta)),a*Mi/M) if (!is.na(Ms)) { if (a*Ms/M < 1) { EB <- rbinom(1,rbinom(1,Ni,(1-px*beta)),a*Ms/M) } else { EB <- rbinom(1,Ni,(1-px*beta)) }} else { EB <- 0 } dNe <- rbinom(1,Ne,omega); dNi <- rbinom(1,Ni,omega); dNr <- rbinom(1,Nr,omega) nNs <- dNe + dNi + dNr Fever <- rbinom(1,Ne,1/lambda) Recover <- rbinom(1,Ni,1/delta) LI <- rbinom(1,Nr,1/eta) Ns <- Ns - IB + nNs + LI Ne <- Ne + IB - dNe - Fever Ni <- Ni - dNi + Fever - Recover Nr <- Nr - dNr + Recover - LI GI <- rbinom(1,Me,1/tau) dMe <- rbinom(1,Me,mu); dMi <- rbinom(1,Mi,mu) nMs <- dMe + dMi Ms <- Ms - EB + nMs; Me <- Me + EB - GI - dMe; Mi <- Mi + GI - dMi } ZI[maxdate] <- Ni/N ZI } execall <- function(px=0.7, beta=0.8, times=50) { PI <- matrix(numeric(maxdate*times), maxdate, times) for (i in 1:times) { PI[,i] <- simexec(px, beta) } matplot(1:maxdate, PI, type="l", col="black", ylim=c(0, 0.4), xlab="days", ylab="parasite rate", main=sprintf("Change of parasite rates as the result of simulation\n [p=%2.1f, beta=%2.1f]", px, beta)) } print(system.time(execall(0.7, 0.8, 50)))