в модели R desolve SIR с петлей и если еще

Question

У меня есть простая модель SIR, я пытаюсь реализовать подход к вакцинации (V), где сначала проверяется, если инфицированные выше порогового значения (100), и если остается достаточно восприимчивых людей (50), он будет вакцинировать определенное количество на временной шаг (50).

Однако я пытаюсь сделать так, чтобы после выполнения условия вакцинация проводилась в течение 7 дней (независимо от того, были ли в течение этих 7 дней инфицированные все еще выше порогового значения или нет, например, если после 4-го дня I = 70). он должен все еще продолжаться, он должен останавливаться только в том случае, если S <50. По истечении 7 дней он должен снова проверить состояние и либо начать снова в течение следующих 7 дней, либо нет.

То, что у меня есть до сих пор, я был бы благодарен, если бы кто-нибудь помог мне реализовать этот цикл

      sirV=function(time, y, params){

S    = y[1]
I    = y[2]
R    = y[3]
V   = y[4]

 with(as.list(params),{
 vac_helper = if (I > 100 & S > 50) {50}

 else {0}
 N = S+I+R+V
 dS = -S*beta*I/N  - vac_helper
 dI = S*beta*I/N - gamma*I
 dR = +gamma*I
 dV = vac_helper

 return(list(c(dS, dI, dR, dV)))
})
}


 myparameters = c(gamma=1/10,beta=0.2)
 times <- seq(0, 300)
 my_ode <- as.data.frame(ode( y=c(100000, 10, 0,0), times, sirV, myparameters))

Answer 1

Вот предложение, но я не совсем уверен, что он ведет себя так, как вы хотели, поэтому сначала проверьте его и вернитесь ко мне, если это не так. Обратите внимание, что end должно быть в глобальной среде.

      library(deSolve)

sirV = function(time, y, params){
  
  S    = y[1]
  I    = y[2]
  R    = y[3]
  V    = y[4]
  
  with(as.list(params),{
    
    # Has the previous vaccination ended, and do we still need to vaccinate?
    if(time > end & I > 100) {
      end <<- time + 7
    }
    
    # Can we vaccinate?
    vaccinate = ifelse(end > time & S > 50, 50, 0)
    
    N = S+I+R+V
    dS = -S*beta*I/N  - vaccinate
    dI = S*beta*I/N - gamma*I
    dR = gamma*I
    dV = vaccinate
    
    # Store some results in a global list so we can check what is happening under the hood
    temp = data.frame(time, end, S, I, R, V, dS, dI, dR, dV, vaccinate)
    catch_results[[length(catch_results)+1]] <<- temp

    return(list(c(dS, dI, dR, dV)))
  })
}

myparameters = c(gamma = 1/10, beta = 0.2)
times <- seq(from = 0, to = 300, by = 1)


end <- 0 # Reset end time
catch_results = list() # Catch results from inside the function
my_ode <- ode( y=c(100000, 10, 0, 0), times, sirV, myparameters)

plot(my_ode)

      # Check this to see if we get the expected behavior, especially at around time = 189
results = dplyr::bind_rows(catch_results)

^{Создано 22.08.2021 пакетом REPEX (v0.3.0)}