NI <- c(3, 8, 28, 75, 221, 291, 255, 235, 190, 125, 70, 28, 12, 5) sir.I <- function(S0=762, I0=1, R0=0, days=1:14, beta=0.0026, gamma=0.5, delta=0.0001) { S <- I <- R <- double(length(days)+1) S[1] <- S0 I[1] <- I0 R[1] <- R0 for (i in days) { S[i+1] <- S[i] - beta*S[i]*I[i] + delta*R[i] I[i+1] <- I[i] + beta*S[i]*I[i] - gamma*I[i] R[i+1] <- R[i] + gamma*I[i] - delta*R[i] } return(I) } plot(0:14, c(1, NI), pch=16, type="p", ylim=c(0, 800), xlab="time (days)", ylab="infected") # model parametes were obtained by data-based estimation # beta comes from the increase of 2 patients on second day # (beta=(I[2]-I[1])/(S[1]*I[1])=2/763=0.0026...), # gamma comes from recovery of most patients after 2 days # (gamma = 1/2). lines(0:14, sir.I(762, 1, 0, 1:14, 0.0026, 0.5, 0.0001), col="red", lwd=2, lty=1) lines(0:14, sir.I(762, 1, 0, 1:14, 0.0013, 0.5, 0.0001), col="green", lwd=2, lty=3) lines(0:14, sir.I(762, 1, 0, 1:14, 0.0039, 0.5, 0.0001), col="black", lwd=2, lty=4) # getting best-fit parameters beta and gamma by Nelder-Mead method params <- c(0.0026, 0.5) diffSIR <- function(z) { mNI <- sir.I(762, 1, 0, 1:14, z[1], z[2], 1/365) RES <- sum((c(1, NI)-mNI)^2) RMSE <- sqrt(RES/15) return(RMSE) } rr <- optim(params, diffSIR, gr="NULL", method="Nelder-Mead") print(rr) lines(0:14, sir.I(762, 1, 0, 1:14, rr\$par[1], rr\$par[2], 1/365), col="blue", lwd=2, lty=2) legend("topright", lty=1:4, col=c("red","blue","green","black"), lwd=2, legend=c("3 cases in day2", "full fit", "2 cases in day2", "4 cases in day2"))