abcdefghij.onestep <- function (x, params) {
Susceptible <- x[2]
Exposed<- x[3]
Infected_Multibacillary <- x[4]
Infected_Paucibacillary<- x[5]
Exposed_Detected_Diagnosis <- x[6]
Treated <- x[7]
Disability<- x[8]
Recovered<- x[9]
Relapse_Multibacillary <-x[10]
Relapse_Paucibacillary<-x[11]
N <- Susceptible + Exposed + Infected_Multibacillary + Infected_Paucibacillary + Exposed_Detected_Diagnosis + Treated + Disability + Recovered + Relapse_Multibacillary + Relapse_Paucibacillary
m12 <- params["m12"]
m25 <- params["m25"]
m23 <- params["m23"]
m24 <- params["m24"]
m35 <- params["m35"]
m45 <- params["m45"]
m37 <- params["m37"]
m56 <- params["m56"]
m67 <- params["m67"]
m68 <- params["m68"]
m89 <- params["m89"]
m810 <- params["m810"]
m96 <- params["m96"]
m97 <- params["m97"]
m106 <- params["m106"]
mu <- params["mu"]
rates <- c(birth=mu*N, susceptible_exposed=m12*Susceptible*Infected_Multibacillary -m25*Infected_Multibacillary * Exposed_Detected_Diagnosis -m23*Exposed*Infected_Multibacillary -m24*Exposed*Infected_Paucibacillary,
exposed_infected_multibacillary=m23*Exposed*Infected_Paucibacillary-m35*Infected_Multibacillary*Exposed_Detected_Diagnosis-m45* Exposed_Detected_Diagnosis* Infected_Paucibacillary-m37*Disability*Infected_Multibacillary,
exposed_infected_paucibacillary=m24*Exposed*Infected_Paucibacillary-m45*Infected_Paucibacillary*Exposed_Detected_Diagnosis,
infected_multibacillary_exposed_detected=m35*Infected_Multibacillary*Exposed_Detected_Diagnosis+m45*Infected_Paucibacillary* Exposed_Detected_Diagnosis-m56* Exposed_Detected_Diagnosis*Treated,
exposed_detected_treatment=m56*Exposed_Detected_Diagnosis*Treated-m67*Treated* Disability-m68*Treated*Recovered,
infected_multibacillary_disability=m37*Disability*Infected_Multibacillary+m67*Treated* Disability,
treated_recovered=m68* Treated*Recovered-m89*Recovered*Infected_Multibacillary-m810* Recovered*Infected_Paucibacillary,
relapse_multibacillary_treatment=-m96*Relapse_Multibacillary* Treated-m97*Relapse_Multibacillary*Disability,
relapse_paucibacillary_treatment=-m106*Relapse_Paucibacillary*Treated,
susceptible_death=mu*Susceptible,
exposed_death=mu*Exposed,
infected_multibacillary_death=mu*Infected_Multibacillary,
infected_paucibacillary_death=mu*Infected_Paucibacillary,
exposed_detected_death=mu* Exposed_Detected_Diagnosis,
treatment_death=mu*Treated,
disability_death=mu*Disability,
recovered_death=mu*Recovered,
relapse_multibacillary_death=mu*Relapse_Multibacillary,
relapse_paucibacillary_death=mu*Relapse_Paucibacillary
)
transitions <- list(birth=c(1,0,0,0,0,0,0,0,0,0),
susceptible_exposed=c(-1,1,0,0,0,0,0,0,0,0),
exposed_infected_multibacillary=c(0,-1,1,0,0,0,0,0,0,0),
exposed_infected_paucibacillary=c(0,-1,0,1,0,0,0,0,0,0),
infected_multibacillary_exposed_detected=c(0,0,-1,0,1,0,0,0,0,0),
exposed_detected_treatment=c(0,0,0,0,-1,1,0,0,0,0),
infected_multibacillary_disability=c(0,0,-1,0,0,0,1,0,0,0),
treated_recovered=c(0,0,0,0,0,-1,0,1,0,0),
relapse_multibacillary_treatment=c(0,0,0,0,0,1,0,0,-1,0),
relapse_paucibacillary_treatment=c(0,0,0,0,0,1,0,0,0,-1),
susceptible_death=c(-1,0,0,0,0,0,0,0,0,0),
exposed_death= c(0,-1,0,0,0,0,0,0,0,0),
infected_multibacillary_death= c(0,0,-1,0,0,0,0,0,0,0),
infected_paucibacillary_death= c(0,0,0,-1,0,0,0,0,0,0),
exposed_detected_death= c(0,0,0,0,-1,0,0,0,0,0),
treatment_death= c(0,0,0,0,0,-1,0,0,0,0),
disability_death= c(0,0,0,0,0,0,-1,0,0,0),
recovered_death= c(0,0,0,0,0,0,0,-1,0,0),
relapse_multibacillary_death= c(0,0,0,0,0,0,0,0,-1,0),
relapse_paucibacillary_death= c(0,0,0,0,0,0,0,0,0,-1)
)
total.rate <- sum(rates)
if (total.rate==0)
tau <- Inf
else
tau <- rexp(n=1,rate=total.rate)
event <- sample.int(n=6,size=1,prob=rates/total.rate)
x+c(tau,transitions[[event]])
}
abcdefghij.simul function
abcdefghij.simul <- function (x, params, maxstep = 10000) {
output <- array(dim=c(maxstep+1,11))
colnames(output) <- names(x)
output[1,] <- x
k <- 1
while ((k <= maxstep) && (x["Exposed"] > 0)) {
k <- k+1
output[k,] <- x <- abcdefghij.onestep(x,params)
}
as.data.frame(output[1:k,])
}
main program
set.seed(56856583)
nsims <- 1
xstart <- c(time=1,Susceptible=100000,Exposed=1,Infected_Multibacillary=1,Infected_Paucibacillary=1,Exposed_Detected_Diagnosis=1,Treated=1,Disability=1,Recovered=1,Relapse_Multibacillary=1,Relapse_Paucibacillary=1)
params <- c(mu=0.02,m12=1,m25=1,m23=1,m24=1,m35=1,m45=1,m37=1,m57=1,m67=1,m68=1,m89=1,m810=1,m96=1,m97=1,m106=1)
library(plyr)
simdat <- rdply(nsims, abcdefghij.simul(xstart,params))
error
Error in if (total.rate == 0) tau <- Inf else tau <- rexp(n = 1, rate = total.rate) : missing value where TRUE/FALSE needed