Advertisement
sin_en

math ia

Nov 15th, 2021
1,589
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 2.95 KB | None | 0 0
  1. run_shiny(model = "SIRS (w/out demography)",
  2.           neweqns = mySIRS,
  3.           ics = c(S = 9999, I = 1, R = 0),
  4.           parm0 = c(beta = 5e-5, gamma = 1/7, delta = 0.1),
  5.           parm_names = c("Transmission rate", "Recovery rate", "Loss of immunity"),
  6.           parm_min = c(beta = 1e-5, gamma = 1/21, delta = 1/9999999),
  7.           parm_max = c(beta = 9e-5, gamma = 1 , delta = 1))
  8.  
  9. run_shiny(model = "mySIR")
  10.           neweqns = mySIR,
  11.           ics = c(S = 9999, I = 1, R = 0),
  12.           parm0 = c(Transmission rate = 1/3, Recovery rate = 1/7),
  13.           parm_names = c("Transmission rate", "Recovery rate"),
  14.           parm_min = c(Transmission rate = 1/30, Recovery rat = 1/21),
  15.           parm_max = c(Transmission rate = 0, Recovery rate = 1))
  16.          
  17. mySIR <- function(t, y, parms) {
  18.  
  19.     with(as.list(c(y, parms)),{
  20.  
  21.         # Change in Susceptibles
  22.         dS <- - beta * S * I
  23.  
  24.         # Change in Infecteds
  25.         dI <- beta * S * I - gamma * I
  26.  
  27.         # Change in Recovereds
  28.         dR <- gamma * I
  29.  
  30.     return(list(c(dS, dI, dR)))
  31.     })
  32. }
  33.  
  34. parm0 = c(R0 = 3, Ip = 7)
  35. parm_names = c("R0", "Infectious period")
  36. parm_min = c(R0 = 0, Ip = 1)
  37. parm_max = c(R0 = 20, Ip = 21
  38.  
  39.  
  40. ## Load deSolve package
  41. library(deSolve)
  42.  
  43. ## Create an SIR function
  44. sir <- function(time, state, parameters) {
  45.  
  46.   with(as.list(c(state, parameters)), {
  47.  
  48.     dS <- -beta * S * I
  49.     dI <-  beta * S * I - gamma * I
  50.     dR <-                 gamma * I
  51.  
  52.     return(list(c(dS, dI, dR)))
  53.   })
  54. }
  55.  
  56. ### Set parameters
  57. ## Proportion in each compartment: Susceptible 0.999999, Infected 0.000001, Recovered 0
  58. init       <- c(S = 1-1e-6, I = 1e-6, R = 0.0)
  59. ## beta: infection parameter; gamma: recovery parameter
  60. parameters <- c(beta = 1.4247, gamma = 0.14286)
  61. ## Time frame
  62. times      <- seq(0, 200, by = 1)
  63.  
  64. ## Solve using ode (General Solver for Ordinary Differential Equations)
  65. out <- ode(y = init, times = times, func = mySIR, parms = parameters)
  66. ## change to data frame
  67. out <- as.data.frame(out)
  68. ## Delete time variable
  69. out$time <- NULL
  70. ## Show data
  71. head(out, 10)
  72.  
  73.  
  74.  
  75.  
  76.  
  77.  mySIR <- function(t, y, parms) {
  78. +
  79. +     with(as.list(c(y, parms)),{
  80. +
  81. +         # Change in Susceptibles
  82. +         dS <- - beta * S * I
  83. +
  84. +         # Change in Infecteds
  85. +         dI <- beta * S * I - gamma * I
  86. +
  87. +         # Change in Recovereds
  88. +         dR <- gamma * I
  89. +
  90. +     return(list(c(dS, dI, dR)))
  91. +     })
  92. + }
  93.  init       <- c(S = 0.999999, I = 0.000001, R = 0.0)
  94.  parameters <- c(beta = 1.1, gamma = 0.3)
  95.  times      <- seq(0, 60, by = 1)
  96.  out <- ode(y = init, times = times, func = mySIR, parms = parameters)
  97.  out <- as.data.frame(out)
  98.  out$time <- NULL
  99.  head(out, 10)
  100. #after data output
  101. matplot(x = times, y = out, type = "l",
  102.         xlab = "Time", ylab = "Susceptible and Recovered", main = "SIR Model",
  103.         lwd = 1, lty = 1, bty = "l", col = 2:4)
  104. legend(40, 0.7, c("Susceptible", "Infected", "Removed"), pch = 1, col = 2:4, bty = "n")
  105.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement