N <- 10 DAYS <- 28 NN <- N*N # CNT <- as.integer(N*(N/2)-N/2) NND <- NN*DAYS # 年齢構造は無視,近接性は考慮 dist2 <- function(X1, X2, Y1, Y2) { (X2-X1)^2+(Y2-Y1)^2 } meandist <- dist2(1, N/2, 1, N/2) # ランダム接触の場合の接触期待値割引きの程度 DISCOUNT <- 2 # matplot()で各人の状態表示。XとYは位置。Sで状態を。 # S = 1は非感染(黒丸),2で感染(赤丸),3で治癒免疫(緑丸) SIMS <- 50 S <- matrix(rep(0, SIMS*DAYS), DAYS, SIMS) I <- matrix(rep(0, SIMS*DAYS), DAYS, SIMS) beta <- 0.5/NN*meandist gamma <- 0.5 palette("R3") pdf("geosims10.pdf", width=14, height=8) set.seed(123) for (s in 1:SIMS) { POP <- data.frame(S = rep(1, NND), X = rep(1:N, N*DAYS), Y = rep(sapply(1:N, rep, N), DAYS)) POP$S[1] <- 2 # POP$S[CNT] <- 2 S[1, s] <- NN - 1 I[1, s] <- 1 for (k in 2:DAYS) { for (i in 1:NN) { if (POP$S[i+(k-2)*NN]==1) { for (j in 1:NN) { betax <- beta/dist2(POP$X[i+(k-2)*NN], POP$X[j+(k-2)*NN], POP$Y[i+(k-2)*NN], POP$Y[j+(k-2)*NN]) if (POP$S[j+(k-2)*NN]==2) { if ((betax>0)&(betax<1)) { POP$S[i+(k-1)*NN] <- 1 + rbinom(1, 1, betax) } } } } else { if (POP$S[i+(k-2)*NN]==2) { POP$S[i+(k-1)*NN] <- 2 + rbinom(1, 1, gamma) } else { POP$S[i+(k-1)*NN] <- POP$S[i+(k-2)*NN] } } } S[k, s] <- sum(POP$S[1:NN+(k-1)*NN]==1) I[k, s] <- sum(POP$S[1:NN+(k-1)*NN]==2) } layout(matrix(1:DAYS, DAYS/7, 7, byrow=TRUE)) for (ii in 1:DAYS) { symbols(POP$X[1:NN+(ii-1)*NN], POP$Y[1:NN+(ii-1)*NN], squares=rep(1, NN), bg=POP$S[1:NN+(ii-1)*NN], inches=FALSE, xlab="", ylab="", main=sprintf("Day %d", ii)) } } layout(1) matplot(1:DAYS, S/NN, type="l", lty=1, lwd=1, col="navy", xlab="DAYS", ylab="Proportion Susceptible", main="Remaining susceptibles with contacts inversely proportional to distance") matplot(1:DAYS, I/NN, type="l", lty=1, lwd=1, col=2, xlab="DAYS", ylab="Proportion Infected", main="Epidemic curve with contacts inversely proportional to distance") set.seed(123) for (s in 1:SIMS) { POP <- data.frame(S = rep(1, NND), X = rep(1:N, N*DAYS), Y = rep(sapply(1:N, rep, N), DAYS)) POP$S[1] <- 2 # POP$S[CNT] <- 2 S[1, s] <- NN - 1 I[1, s] <- 1 for (k in 2:DAYS) { for (i in 1:NN) { if (POP$S[i+(k-2)*NN]==1) { for (j in 1:NN) { betax <- beta/DISCOUNT if (POP$S[j+(k-2)*NN]==2) { if ((betax>0)&(betax<1)) { POP$S[i+(k-1)*NN] <- 1 + rbinom(1, 1, betax) } } } } else { if (POP$S[i+(k-2)*NN]==2) { POP$S[i+(k-1)*NN] <- 2 + rbinom(1, 1, gamma) } else { POP$S[i+(k-1)*NN] <- POP$S[i+(k-2)*NN] } } } S[k, s] <- sum(POP$S[1:NN+(k-1)*NN]==1) I[k, s] <- sum(POP$S[1:NN+(k-1)*NN]==2) } layout(matrix(1:DAYS, DAYS/7, 7, byrow=TRUE)) for (ii in 1:DAYS) { symbols(POP$X[1:NN+(ii-1)*NN], POP$Y[1:NN+(ii-1)*NN], squares=rep(1, NN), bg=POP$S[1:NN+(ii-1)*NN], inches=FALSE, xlab="", ylab="", main=sprintf("Day %d", ii)) } } layout(1) matplot(1:DAYS, S/NN, type="l", lty=1, lwd=1, col="navy", xlab="DAYS", ylab="Proportion Susceptible", main="Remaining susceptible with random contacts") matplot(1:DAYS, I/NN, type="l", lty=1, lwd=1, col=2, xlab="DAYS", ylab="Proportion Infected", main="Epidemic curve with random contacts") dev.off()