dajmcdon / rtestim

https://dajmcdon.github.io/rtestim/
Other
4 stars 0 forks source link

Return negative total infectiousness for the first few observations #60

Closed jiapivialiu closed 3 months ago

jiapivialiu commented 4 months ago

Run the following code the reproduce the error:

length <- 300
## piecewise constant Rt: 
Rt <- c(rep(2, floor(2 * length / 5)), rep(0.8, length - floor(2 * length / 5)))
Mu <- double(length); Mu[1] <- 2
incidence <- double(length)
w <- double(length)
## gamma parameters of measles: 
gamma_pars <- c(14.5963, 1.0208) 

set.seed(317)
incidence[1] <- Mu[1]
for(t in 2:length){
  w[t] <- delay_calculator(incidence[1:t], dist_gamma = gamma_pars)[t]
  Mu[t] <- Rt[t] * w[t]
  incidence[t] <- rpois(1, Mu[t])
}
w_full <- delay_calculator(incidence, dist_gamma = gamma_pars)
w_full[1:10]
w[1:10]

The second component of w, i.e., the total infectiousness, is generated to be negative, which causes the bug.

dajmcdon commented 4 months ago

Note: The result is

[1] -0.002660358 -0.002660358  2.000004916  2.013321195  2.073505696
[6]  2.181887784  2.309699224  2.432719834  2.547963903  2.671985185