—
title: “Lab6”
author: “Jiaqi Gao”
date: “11/20/2019”
output: html_document
—
# imput data
“`{r}
times=c(0,7,14,21,28,35,42,49,56,63,70,77,84,91,98)
data=c(97,271,860,1995,4419,6549,6321,4763,2571,1385,615,302,159,72,34)
“`
# 1) SIR model set up
“`{r}
library(deSolve)
SIR_ODE=function(t,x,param){
B=param[1]
y=param[2]
k=param[3]
S=x[1]
I=x[2]
R=x[3]
dxdt=numeric(3)
dxdt[1]=-B*S*I
dxdt[2]=B*S*I-y*I
dxdt[3]=y*I
return(list(dxdt))
}
# times for the simulation to run on
times=seq(0,98,7)# 0 to 100 days, counting by 1 day
# initial values
x0=c(“S”=0.999, “I”=0.001, “R”=0)
# set up parameter values
paramvals=c(“B”=0.4, “y”=0.25, “k”=1)
# simulate the model
xsim1=ode(x0, times, SIR_ODE, paramvals, method=’ode45′)
# formulate Y
N=200000
k=1
Y_k1=k*N*xsim1[,3]
# plot Y
plot(xsim1[,1], Y_k1, type=’l’,xlab=”days”,ylab=”number of susceptible people get infected”,col=”brown”)
points(x=xsim1[,1], y=data, type=”p”)
“`
# calculate likelihood
“`{r}
Cal_SIR=function(param,SIR_ODE, x0, times, data){
param=abs(param)
xcurr=ode(x0, times,SIR_ODE, param)
k=param[3]
Y_k1=k*N*xcurr[,3]
N=200000
NLL_SIR=sum(Y_k1)-sum(data*log(Y_k1))
return(NLL_SIR)
}
times=c(0,7,14,21,28,35,42,49,56,63,70,77,84,91,98)
x0=c(“S”=0.999, “I”=0.001, “R”=0)
data=c(97,271,860,1995,4419,6549,6321,4763,2571,1385,615,302,159,72,34)
param1=c(“B”=0.4, “y”=0.25, “k”=1)
param2=c(“B”=0.38, “y”=0.25, “k”=1)
param3=c(“B”=0.42, “y”=0.25, “k”=1)
Cal_SIR(param1,SIR_ODE, x0, times, data)
Cal_SIR(param2,SIR_ODE, x0, times, data)
Cal_SIR(param3,SIR_ODE, x0, times, data)
“`
# 3) parameter estimates
“`{r}
x0=c(“S”=0.999, “I”=0.001, “R”=0)
times=c(0,7,14,21,28,35,42,49,56,63,70,77,84,91,98)
data=c(97,271,860,1995,4419,6549,6321,4763,2571,1385,615,302,159,72,34)
results=optim(param1, fn=Cal_SIR, SIR_ODE=SIR_ODE, times=times, data=data, x0=x0)
parameters=results$par
negLL=results$value
print(results)
“`
# plot the model and the data
“`{r}
SIR_ODE1=function(t,x,param){
B=param[1]
y=param[2]
k=param[3]
S=x[1]
I=x[2]
R=x[3]
dxdt=numeric(3)
dxdt[1]=-B*S*I
dxdt[2]=B*S*I-y*I
dxdt[3]=y*I
return(list(dxdt))
}
# times for the simulation to run on
times=seq(0,98,7)# 0 to 100 days, counting by 1 day
# initial values
x0=c(“S”=0.999, “I”=0.001, “R”=0)
# set up parameter values
paramvals=c(“B”=0.44, “y”=0.29, “k”=0.54)
# simulate the model
xsim2=ode(x0, times, SIR_ODE, paramvals, method=’ode45′)
# formulate Y
N=200000
k=0.54
Y_k0.54=k*N*xsim2[,3]
# plot Y
plot(xsim2[,1], Y_k0.54, type=’l’,xlab=”days”,ylab=”number of susceptible people get infected”,col=”brown”)
points(x=xsim2[,1], y=data, type=”p”)
“`