程序代写代做代考 —


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”)
“`