# http://minato.sip21c.org/tiid/malsim.R 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 } PI <- simexec(0.7,0.8) plot(1:maxdate,PI,type="l",ylim=c(0,0.4),xlab="days",ylab="parasite rate", main="Change of parasite rates as the result of simulation\n [p=0.7, beta=0.8]") for (i in 2:50) { PI <- simexec(0.7,0.8) lines(1:maxdate,PI) }