noah/CS-7863-Sci-Stat-Final
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Model v1, with call to python

Browse Source
This commit is contained in:
Noah L. Schrick 2023-05-01 15:03:31 -05:00
parent 7ac37f0aea
commit 184be47a33
3 changed files with 27 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
Analysis/SERIDS.pdf Normal file
View File

Binary file not shown.

50
Analysis/ode.R
View File

@ -17,11 +17,11 @@ seirds.f <- function(t, y, k) {
waning = k[7] waning = k[7]
# Saying infec rate of S in contact with E same as contact with I # Saying infec rate of S in contact with E same as contact with I
dS <- epsilon - (beta*E + beta*I)*S + waning*R - mu*S dS <- epsilon - (beta*E + beta*I)*S + waning*R - gamma_d*S
dE <- beta*S*E + beta*S*I - (delta+mu)*E dE <- beta*S*E + beta*S*I - (delta+gamma_d)*E
dI <- delta*E - mu*I - gamma_r*I dI <- delta*E - (1+gamma_d)*I
dR <- gamma_r*I - (waning+mu)*R dR <- (1-mu)*I - (waning+gamma_d)*R
dD <- gamma_d*I dD <- mu*I
return(as.matrix(c(dS, dE, dI, dR, dD))) return(as.matrix(c(dS, dE, dI, dR, dD)))
} }
@ -33,30 +33,23 @@ library(reticulate)
#conda_install(envname = "r-reticulate", packages="pygraphviz") #conda_install(envname = "r-reticulate", packages="pygraphviz")
source_python('prep_model.py') source_python('prep_model.py')
# a) a=0.5, b=1, S(0)=0.9, I(0)=0.1, R(0)=0 # Obtained from prep_model.py
seirds.params <- c(0.7316455696202532, # beta seirds.params <- c(beta, # beta
0.020253164556962026, # delta delta, # delta
0.09113924050632911, # gamma_r gamma_r, # gamma_r
0.002531645569620253, # gamma_d gamma_d, # gamma_d
0.0034602076124567475, # mu mu, # mu
0.0, # epsilon epsilon, # epsilon
0.0) # waning omega) # waning
S <- 62 S <- S/length(color_d)
E <- 8 E <- E/length(color_d)
I <- 288 I <- (I_R+I_D)/length(color_d)
R <- 36 R <- R/length(color_d)
D <- 1 D <- D/length(color_d)
N <- S+E+I+R+D
S <- S/N
E <- E/N
I <- I/N
R <- R/N
D <- D/N
tmin <- 0 tmin <- 0
tmax <- 7 tmax <- 5
library(pracma) library(pracma)
seirds.ode.sol <- ode45(f = function(t,y){seirds.f(t,y,k=seirds.params)}, seirds.ode.sol <- ode45(f = function(t,y){seirds.f(t,y,k=seirds.params)},
@ -86,5 +79,6 @@ plot.seirds <- function(sol, method){
} }
plot.seirds(seirds.ode.sol, "ODE45") plot.seirds(seirds.ode.sol, "ODE45")
ggsave("SERIDS.pdf")
print(seirds.ode.sol$y[tmax,1]+seirds.ode.sol$y[tmax,2]+seirds.ode.sol$y[tmax,3]+seirds.ode.sol$y[tmax,4]+seirds.ode.sol$y[tmax,5]) # Sanity check: Make sure sums to ~1.0
sum(tail(seirds.ode.sol$y,1))

7
Analysis/prep_model.py
View File

@ -85,12 +85,13 @@ for node in A:
# Params # Params
beta = (I_R+I_D)/len(A) # rate of infec (I/total?) beta = (I_R+I_D)/len(A) # rate of infec (I/total?)
delta = E/len(A) # symptom appearance rate (E/total?) #delta = E/len(A) # symptom appearance rate (E/total?)
delta = 1 # incubation period
gamma_r = R/len(A) # recov rate (R/total?) gamma_r = R/len(A) # recov rate (R/total?)
gamma_d = D/len(A) # death rate (D/total?) gamma_d = D/len(A) # death rate (D/total?)
mu = D/(I_R+I_D) # fatality ratio (D/I) mu = D/(I_R+I_D) # fatality ratio (D/I)
epsilon = ep_tmp/len(A) # infected import rate epsilon = ep_tmp/len(A) # infected import rate
omega = 0 # waning immunity rate omega = 1 # waning immunity rate
print("Model Compartments:") print("Model Compartments:")
@ -109,3 +110,5 @@ print("gamma_r:", str(gamma_r))
print("gamma_d:", str(gamma_d)) print("gamma_d:", str(gamma_d))
print("mu:", str(mu)) print("mu:", str(mu))
print("epsilon:", str(epsilon)) print("epsilon:", str(epsilon))
print("omega:", str(omega))