-
Notifications
You must be signed in to change notification settings - Fork 0
/
cardsim.R
29 lines (27 loc) · 1.07 KB
/
cardsim.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# Transmission function, for a given number of contacts in a particular model world. det=TRUE returns the deterministic version; cards=TRUE a stochastic version based on the detailed model; cards=FALSE a binomial approximation to the "true" stochastic version. They should be very similar, and produce qualitatively similar results
transFun <- function(contacts, S, N, det=TRUE, cards=TRUE){
if(contacts==0) return(0)
if(contacts>=N) return(S)
if (det) return(contacts*S/N)
## Pick contacts from population, and see how many are in the first S
if (!cards) return(rbinom(1, size=S, prob=contacts/N))
return(sum(sample(N, contacts)<=S))
}
# Simulate a stochastic epidemic, using the transFun above for the transmission step
sim <- function(R0, N=26, numSteps=10, I_0=1, det=TRUE, cards=TRUE){
Ivec <- I <- I_0
Svec <- S <- N-I_0
for (j in 2:numSteps){
trans <- transFun(R0*I, S, N, det, cards)
I <- trans
S <- S-trans
if (S<0) print(c(R0*I, S, N, det, cards, trans))
Ivec[[j]] <- I
Svec[[j]] <- S
}
return(data.frame(
timestep=1:numSteps
, I=Ivec
, S=Svec
))
}