diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5b6a065 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +.Rproj.user +.Rhistory +.RData +.Ruserdata diff --git a/Analysis/ode.R b/Analysis/ode.R new file mode 100644 index 0000000..dd00951 --- /dev/null +++ b/Analysis/ode.R @@ -0,0 +1,87 @@ +seirds.f <- function(t, y, k) { + S <- y[1] + E <- y[2] + I_R <- y[3] + I_D <- y[4] + R <- y[5] + D <- y[6] + N <- y[7] + + beta = k[1] + delta = k[2] + gamma_r = k[3] + gamma_d = k[4] + mu = k[5] + epsilon = k[6] + waning = k[7] + + dS <- S - beta*(I_R+I_D)/N + waning*R + dE <- E + beta*S*(I_R+I_D)/N - delta*E + epsilon + dI_R <- I_R + delta*(1-mu)*E - gamma_r*I_R + epsilon + dI_D <- I_D + delta*mu*E - gamma_d*I_D + epsilon + dR <- R + gamma_r*I_R - waning*R + dD <- D + gamma_d*I_D + return(as.matrix(c(dS, dE, dI_R, dI_D, dR, dD, N))) +} + +# a) a=0.5, b=1, S(0)=0.9, I(0)=0.1, R(0)=0 +seirds.params <- c(0.7316455696202532, # beta + 0.020253164556962026, # delta + 0.09113924050632911, # gamma_r + 0.002531645569620253, # gamma_d + 0.0034602076124567475, # mu + 0.0) #epsilon + +seirds.params <- c(1, # beta + 1, # delta + 1, # gamma_r + 1, # gamma_d + 1, # mu + 1) #epsilon + + +S <- 62 +E <- 8 +I_R <- 288 +I_D <- 1 +R <- 36 +D <- 1 + +N <- S+E+I_R+I_D+R+D +S <- S/N +E <- E/N +I_R <- I_R/N +I_D <- I_D/N +R <- R/N +D <- D/N + +tmin <- 0 +tmax <- 20 + +install.packages("pracma") +library(pracma) +seirds.ode.sol <- ode45(f = function(t,y){seirds.f(t,y,k=seirds.params)}, + y = c(S, E, I_R, I_D, R, D, N), + t0 = tmin, tfinal = tmax) +# plot +plot.seirds <- function(sol, method){ + # - pred/prey columns melted to 1 column called "value" + sol.melted <- melt(data.frame(time=sol$t, + S=sol$y[,1], + E=sol$y[,2], + I_R=sol$y[,3], + I_D=sol$y[,4], + R=sol$y[,5], + D=sol$y[,6],), + id = "time") # melt based on time + # New column created with variable names, called “variable” + colnames(sol.melted)[2] <- "Group" # used in legend + g <- ggplot(data = sol.melted, + aes(x = time, y = value, color = Group)) + g <- g + geom_point(size=2) + g <- g + xlab("time") + ylab("Population") + g <- g + ggtitle(paste("SEIRDS Model Using", method)) + show(g) +} + +plot.seirds(seirds.ode.sol, "ODE45") \ No newline at end of file diff --git a/CS-7863-Sci-Stat-Final.Rproj b/CS-7863-Sci-Stat-Final.Rproj new file mode 100644 index 0000000..8e3c2eb --- /dev/null +++ b/CS-7863-Sci-Stat-Final.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX