炎上案件を見かけた。本当は祭りになっているときに早めに投下して鮮度を大事にしたかったけど、個人的用事で忙しいこととシミュレーションプログラムを作るのに時間がかかったこともあり、この記事の鮮度はもう0よｯ…
おそらく一般人が全然知らないワクチンの効能として、ワクチンを打った自分自身が病気にかかりにくくなる、というのは当然のこととして、他人が病気にかかりにくくなり、集団全体として大流行を防ぐことができる、というのがある。これをSIR モデルという数理モデルで情報科学的にR によるシミュレーションで示そう。
教室に規則正しく並べた机のごとく、の大きさの適当な正方形(長方形でもいいけど)を考える。1セルが1人の人に相当するし、正方形全体がとある集団、村とか都市とか、そういうのに相当する。
ここに、外から流入してきたとか、自然発生したとか、時刻0 の時点で何人かが感染症を発症したとする(start_infect)。この感染症は、とある感染者1人の周囲8セルに確率論的に[tex:p\hspace{3}(0
nr <- 150 # 自乗が人口 i_day <- 3 # 感染力のある日数 obs_time <- 100 # 観察期間 start_infect <- 20 # 最初の感染者数 infect_prob <- 0.7 # 感染力 vaccine_power <- 0.9 # ワクチンによって感染力がどれくらい下がるかの効果 vaccine_prop <- 0.950 # ワクチン接種している人の割合 mat <- mat0 <- v <- diag(0, nr+2) mat[c(1, nr+2),] <- 1 mat[,c(1, nr+2)] <- 1 idx0 <- which(c(mat) == 1) # 周辺のいらない項 idx1 <- seq((nr+2)^2)[-idx0] # 実際の人の項 mat[idx0] <- 0 vaccine_n <- floor(vaccine_prop*nr^2) # ワクチン接種している人数 v[sample(idx1, vaccine_n)] <- 1 # ワクチン接種した人 ### sir <- list(s=v, i=mat, r=mat, v=mat) idxr <- which(sir$r > 0) # 既感染の id sir$i[sample(setdiff(idx1, which(v == 1)), start_infect)] <- 1 # 最初の感染者 res <- matrix(0, nr=obs_time, nc=length(sir)) colnames(res) <- names(sir) b <- matrix(0, nr=obs_time, nc=5) xl <- c(1, obs_time) yl <- c(0, 1) library(animation) #saveGIF( for(j in seq(obs_time)){ infected <- setdiff(which(sir$i > 0), idxr) # 現感染者の id idx8 <- unlist(cbind(mapply(function(x) infected+x, -1:1), cbind.data.frame(mapply(function(k) mapply(function(j) infected+(nr+j)*k, 1:3), c(1, -1), SIMPLIFY=FALSE))), use.name=FALSE) # 現感染者の周囲8人 idx8 <- idx8[!(idx8 %in% union(idx0, c(infected, idxr)))] # 周辺項と既感染者を除く a <- ifelse( tapply(mapply(function(x) sample(0:1, size=1, prob=c(1, 0)+c(-1, 1)*infect_prob*ifelse(x==1, 1-vaccine_power, 1)), v[idx8]), idx8, sum)>0, 1, 0) # 現感染者から周囲に感染したか sir$i <- ifelse(sir$i>0, sir$i+1, 0) # 感染している人に感染時間を足す sir$i[as.numeric(names(a))[a>0]] <- 1 # 新たに感染した人 sir$r[sir$i > i_day] <- 1 # 既感染 idxr <- which(sir$r > 0) # 既感染の id sir$v[which((sir$i>0) * (v>0) > 0)] <- 1 # ワクチン打ったのに感染した人 k <- mapply(function(x) ifelse(x[idx1]>1, 1, x[idx1]), sir) # 人数の計算 k[,"i"] <- replace(k[,"i"], k[,"r"]==1, 0) # 現感染者と既感染者をかぶらないように res[j,] <- colMeans(k) b[j,] <- c(sum(res[j,c("i","r")]), res[j,c("i","r")], mean(k[k[,"s"]==1,"v"]), mean(k[k[,"s"]==0, "r"])) flag <- rowSums(sweep(k[,-4] ,2, 2^(0:2), "*")) cols <- c("white", "lightgreen", "red", "yellow", "black", grey(0.8)) # ワクチンを打たずに感染してない人 # ワクチンを打ってて感染してない人 # ワクチンを打たずに感染した人 # ワクチンを打ってて感染した人 # ワクチンを打たずに既感染の人 # ワクチンを打ってて既感染の人 mat0[idx1] <- flag mat0[1, 1:length(cols)] <- seq(cols)-1 xy <- seq(nr+2) par(mfrow=c(1, 2)) image(xy, xy, mat0, col=cols, frame=F, axes=F, xlab="", ylab="") title(paste("Vaccination", round(vaccine_prop, 3))) polygon(c(0.5,1.5,1.5,0.5), rep(c(0.5, length(cols)+0.5), each=2), col="white", border=NA) matplot(b[seq(j),], type="l", las=1, lwd=3, xlim=xl, ylim=yl, lty=1, xlab="Time", ylab="Proportion") legend("topleft", legend=c("総感染者", "現感染者", "既感染者", "ワクチン感染者", "非ワクチン感染者"), lwd=3, col=seq(ncol(b))) } #, movie.name = "brownian_motion.gif", interval = 0.1, nmax = 30, ani.width = 800, ani.height = 400)
