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) } windowsFonts(JP1=windowsFont("MS Gothic"), JP2=windowsFont("MS Mincho"), JP3=windowsFont("Meiryo")) par(family="JP3") 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 # gamma comes from recovery of most patients after 2 days lines(0:14, sir.I(762, 1, 0, 1:14, 0.0026, 0.5, 0.0001), col="red", lty=1) lines(0:14, sir.I(762, 1, 0, 1:14, 0.0013, 0.5, 0.0001), col="green", lty=1) lines(0:14, sir.I(762, 1, 0, 1:14, 0.0039, 0.5, 0.0001), col="darkgray", lty=1) # 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", lty=2) legend("top", lty=c(1,1,1,2), col=c("red","green","darkgray","blue"), legend=c("2日目の発症2人（β=0.0026）","もし1人なら（β=0.0013）","もし3人なら（β=0.0039）","全データへの当てはめ"))