Final project

set.seed(1345)

Tmax <- 100
N <- 1000

p <- 0.0015
q <- 1 - p

S <- numeric(Tmax + 1)
I <- numeric(Tmax + 1)
A <- numeric(Tmax + 1)

S[1] <- 998
I[1] <- 2
A[1] <- 0

for (t in 1:Tmax) {
  prob_infect <- 1 - q^I[t]
  new_inf <- rbinom(1, S[t], prob_infect)

  I[t + 1] <- new_inf
  S[t + 1] <- S[t] - new_inf
  A[t + 1] <- A[t] + new_inf
}

barplot(I,
        xlab = "Time",
        ylab = "New infections",
        main = "Reed Frost Simulation R0 = 1.5")