R语言模拟疫情传播-RVirusBroadcast
原创：hxj7
本文用RVirusBroadcast展示模拟的疫情数据
本文篇幅较长，分为以下几个部分：
- 前言
- 效果展示
- 小结
- 附录：RVirusBroadcast代码
前言
前几天微博的一个热搜主题是“计算机仿真程序告诉你为什么现在还没到出门的时候！！！”，该视频用模拟的疫情数据告诉大家“不要随便出门（宅在家）”对战胜疫情很重要，生动形象，广受好评。
所用的程序叫VirusBroadcast，源码已经公开，是用Java写的。鉴于画图是R语言的优势，所以笔者在读过源码后，写了一个VirusBroadcast程序的R语言版本，暂且叫做RVirusBroadcast。与VirusBroadcast相比，RVirusBroadcast所用的模型和逻辑大体不变，只是在少许细节上做了修改。
（为了防止上面的超链接被过滤掉而打不开，文末也放上了明文链接）
效果展示
下面两段视频是RVirusBroadcast用模拟的数据展示的效果，由于笔者的电脑性能实在一般，所以暂时只模拟了30天的数据。请再次注意下面两段视频的数据是模拟生成的，纯属虚构，不具有现实意义，仅供电脑模拟实验所用。
其他条件不变，当人们随意移动时，病毒传播迅速，疫情很难控制
R语言模拟疫情传播（随意移动）https://www.zhihu.com/video/1211607442850734080其他条件不变，当人们控制自己的移动时，病毒传播缓慢，疫情逐渐得到控制
R语言模拟疫情传播（控制移动）https://www.zhihu.com/video/1211607589467205632小结
诚如VirusBroadcast的作者所说，现在的模型是一个很简单的模型，所用的数据也是模拟生成的，还需优化改进。朋友们如果有兴趣，可以自行查阅复制下文中的R代码，自由修改。
如果您对代码有任何意见或建议，请联系hxj5hxj5@126.com。谢谢！
参考
[1] “计算机仿真程序告诉你为什么现在还没到出门的时候” 视频地址：
https://www.bilibili.com/video/av86478875?spm_id_from=333.788.b_765f64657363.1
[2] VirusBroadcast (Java)程序源码：
https://github.com/KikiLetGo/VirusBroadcast
附录：RVirusBroadcast代码
1 ###name:RVirusBroadcast
2 ###author:hxj7(hxj5hxj5@126.com)
3 ###version:202002010
4 ###note:本程序是"VirusBroadcast (in Java)"的R版本
5 ### VirusBroadcast (in Java) 项目链接：
6 ### https://github.com/KikiLetGo/VirusBroadcast/tree/master/src
7
8 library(tibble)
9 library(dplyr)
10
11 ########## 模拟参数 ##########
12 ORIGINAL_COUNT <- 50 # 初始感染数量
13 BROAD_RATE <- 0.8 # 传播率
14 SHADOW_TIME <- 140 # 潜伏时间，14天为140
15 HOSPITAL_RECEIVE_TIME <- 10 # 医院收治响应时间
16 BED_COUNT <- 1000 # 医院床位
17
18 MOVE_WISH_MU <- -0.99 # 流动意向平均值，建议调整范围：[-0.99,0.99];
19 # -0.99 人群流动慢速率，甚至完全控制疫情传播;
20 # 0.99为人群流动快速率, 可导致全城感染
21
22 CITY_PERSON_SIZE <- 5000 # 城市总人口数量
23
FATALITY_RATE <- 0.02 # 病死率，根据2月6日数据估算（病死数/确诊数）为0.02
24 SHADOW_TIME_SIGMA <- 25 # 潜伏时间方差
25 CURED_TIME <- 50 # 治愈时间均值，从入院开始计时
26 CURED_SIGMA <- 10 # 治愈时间标准差
27 DIE_TIME <- 300 # 死亡时间均值，30天，从发病（确诊）时开始计时
28 DIE_SIGMA <- 50 # 死亡时间标准差
29
30 CITY_WIDTH <- 700 # 城市大小即窗口边界，限制不允许出城
31 CITY_HEIGHT <- 800
32
33 MAX_TRY <- 300 # 大模拟次数，300代表30天
34
35 ########## 生成人群点，用不同颜色代表不同健康状态。 ##########
36 # 用正态分布刻画人群点的分布
37 CITY_CENTERX <- 400 # x轴的mu值
38 CITY_CENTERY <- 400
39 PERSON_DIST_X_SIGMA <- 100 # x轴的sigma值
40 PERSON_DIST_Y_SIGMA <- 100
41
42 # 市民状态应该需要细分，虽然有的状态暂未纳入模拟，但是细分状态应该保留
43 STATE_NORMAL <- 0 # 正常人，未感染的健康人
44 STATE_SUSPECTED <- STATE_NORMAL + 1 # 有暴露感染风险
45 STATE_SHADOW <- STATE_SUSPECTED + 1 # 潜伏期
46 STATE_CONFIRMED <- STATE_SHADOW + 1 # 发病且已确诊为感染病人
47 STATE_FREEZE <- STATE_CONFIRMED + 1 # 隔离治疗，禁止位移
48 STATE_DEATH <- STATE_FREEZE + 1 # 病死者
49 STATE_CURED <- STATE_DEATH + 1 # 治愈数量用于计算治愈出院后归还床位数量，该状态是否存续待定
50
51 worldtime <- 0
52 NTRY_PER_DAY <- 10 # 一天模拟几次
53 getday <- function(t) (t - 1) %/% NTRY_PER_DAY + 1
54
55 # 生成人群数据
56 format_coord <- function(coord, boundary) {
57 if (coord < 0) return(runif(1, 0, 10))
58 else if (coord > boundary) return(runif(1, boundary - 10, boundary))
59 else return(coord)
60 }
61 set.seed(123)
62 people <- tibble(
63 id = 1:CITY_PERSON_SIZE,
64 x = sapply(rnorm(CITY_PERSON_SIZE, CITY_CENTERX, PERSON_DIST_X_SIGMA),
65 format_coord, boundary = CITY_WIDTH), # (x, y) 为人群点坐标
66 y = sapply(rnorm(CITY_PERSON_SIZE, CITY_CENTERY, PERSON_DIST_Y_SIGMA),
67 format_coord, boundary = CITY_HEIGHT),
68 state = STATE_NORMAL, # 健康状态
69 infected_time = 0, # 感染时刻
70 confirmed_time = 0, # 确诊时刻
71 freeze_time = 0, # 隔离时刻
72 cured_moment = 0, # 痊愈时刻，为0代表不确定
73 die_moment = 0 # 死亡时刻，为0代表未确定，-1代表不会病死
74 ) %>%
75 mutate(tx = rnorm(CITY_PERSON_SIZE, x, PERSON_DIST_X_SIGMA), # target x
76 ty = rnorm(CITY_PERSON_SIZE, y, PERSON_DIST_Y_SIGMA),
77 has_target = T, is_arrived = F)
78
79# 随机选择初始感染者
80 peop_id <- sample(people$id, ORIGINAL_COUNT)
81 people$state[peop_id] <- STATE_SHADOW
82 people$infected_time[peop_id] <- worldtime
83 people$confirmed_time[peop_id] <- worldtime +
84 max(rnorm(length(peop_id), SHADOW_TIME / 2, SHADOW_TIME_SIGMA), 0)
85
86 ########## 生成床位点 ##########
87 HOSPITAL_X <- 720 # 张床位的x坐标
88 HOSPITAL_Y <- 80 # 张床位的y坐标
89 NBED_PER_COLUMN <- 100 # 医院每一列有多少张床位
90 BED_ROW_SPACE <- 6 # 一行中床位的间距
91 BED_COLUMN_SPACE <- 6 # 一列中床位的间距
92
93 bed_ncolumn <- ceiling(BED_COUNT / NBED_PER_COLUMN)
94 hosp_beds <- tibble(id = 1, x = 0, y = 0, is_empty = T, state = STATE_NORMAL) %>%
95 slice(-1)
96 if (BED_COUNT > 0) {
97 hosp_beds <- tibble(
98 id = 1:BED_COUNT,
99 x = HOSPITAL_X + rep(((1:bed_ncolumn) - 1) * BED_ROW_SPACE,
100 each = NBED_PER_COLUMN)[1:BED_COUNT],
101 y = HOSPITAL_Y + 10 - BED_COLUMN_SPACE +
102 rep((1:NBED_PER_COLUMN) * BED_COLUMN_SPACE, bed_ncolumn)[1:BED_COUNT],
103 is_empty = T,
104 person_id = 0 # 占用床位的患者的序号，床位为空时为0
105 )
106 }
107
108 ########## 准备画图的数据 ##########
109 npeople_total <- CITY_PERSON_SIZE
110 npeople_shadow <- ORIGINAL_COUNT
111 npeople_confirmed <- npeople_freeze <- npeople_cured <- npeople_death <- 0
112 nbed_need <- 0
113
114 ########## 画出初始数据 ##########
115 # 设置画图参数
116 person_color <- data.frame( # 不同健康状态的颜色不同
117 label = c("健康", "潜伏", "确诊", "隔离", "治愈", "死亡"),
118 state = c(STATE_NORMAL, STATE_SHADOW, STATE_CONFIRMED, STATE_FREEZE,
119 STATE_CURED, STATE_DEATH),
120 color = c(
121 "lightgreen", # 健康
122 "#EEEE00", # 潜伏期
123 "red", # 确诊
124 "#FFC0CB", # 隔离
125 "green", # 治愈
126 "black" # 死亡
127 ), stringsAsFactors = F
128 )
129 bed_color <- data.frame(
130 is_empty = c(T, F), color = c("#F8F8FF", "#FFC0CB"), stringsAsFactors = F
131 )
132 x11(width = 5, height = 7, xpos = 0, ypos = 0, title = "人群变化模拟")
133 window_hist <- dev.cur()
134 x11(width = 7, height = 7, xpos = 460, ypos = 0, title = "疫情传播模拟")
135 window_scatter <- dev.cur()
136 max_plot_x <- ifelse(BED_COUNT > 0, max(hosp_beds$x), CITY_WIDTH) + 10
137
138 # 疫情传播模拟散点图
139 dev.set(window_scatter)
140 plot(x = people$x, y = people$y, cex = .8, pch = 20, xlab = NA, ylab = NA,
141 xlim = c(5, max_plot_x), xaxt = "n", yaxt = "n", bty = "n", main = "疫情传播模拟",
142 sub = paste0("世界时间第 ", getday(worldtime), " 天"),
143 col = (people %>% left_join(person_color, by = "state") %>%
144 select(color))$color)
145 points(x = hosp_beds$x, y = hosp_beds$y, cex = .8, pch = 20,
146 col = (hosp_beds %>% left_join(bed_color, by = "is_empty") %>%
147 select(color))$color)
148 rect(HOSPITAL_X - BED_ROW_SPACE / 2, HOSPITAL_Y + 10 - BED_COLUMN_SPACE,
149 max(hosp_beds$x) + BED_ROW_SPACE / 2, max(hosp_beds$y + BED_COLUMN_SPACE))
150 legend(x = 150, y = -30, legend = person_color$label, col = person_color$color,
151 pch = 20, horiz = T, bty = "n", xpd = T)
152
153 # 人群变化模拟条形图
154 dev.set(window_hist)
155 bp_data <- c(npeople_death, npeople_cured, nbed_need, npeople_freeze,
156 npeople_confirmed, npeople_shadow)
157 bp_color <- c("black", "green", "#FFE4E1", "#FFC0CB", "red", "#EEEE00")
158 bp_labels <- c("死亡", "治愈", "不足\n床位", "隔离", "累计\n确诊", "潜伏")
159 bp <- barplot(bp_data, horiz = T, border = NA, col = bp_color,
160 xlim = c(0, CITY_PERSON_SIZE * 1), main = "人群变化模拟",
161 sub = paste0("世界时间第 ", getday(worldtime), " 天"))
162 abline(v = BED_COUNT, col = "gray", lty = 3)
163 abline(v = CITY_PERSON_SIZE, col = "gray", lty = 1)
164 text(x = -350, y = bp, labels = bp_labels, xpd = T)
165 text(x = bp_data + CITY_PERSON_SIZE / 15, y = bp, xpd = T,
166 labels = ifelse(bp_data > 0, bp_data, ""))
167 legend(x = 300, y = -.6, legend = c("总床位数", "城市总人口"), col = "gray",
168 lty = c(3, 1), bty = "n", horiz = T, xpd = T)
169
170 Sys.sleep(5) # 手动调整窗口大小
171
172 ########## 更新人群数据 ##########
173 # 市民流动意愿以及移动位置参数174 MOVE_WISH_SIGMA <- 1
175 MOVE_DIST_SIGMA <- 50
176 SAFE_DIST <- 2 # 安全距离
177
178 worldtime <- worldtime + 1
179 get_min_dist <- function(person, peop) { # 一个人和一群人之间的小距离
180 min(sqrt((person["x"] - peop$x) ^ 2 + (person["y"] - peop$y) ^ 2))
181 }
182 for (i in 1:MAX_TRY) {
183 # 如果已经隔离或者死亡了，就不需要处理了
184 #
185 # 处理已经确诊的感染者（即患者）
186 peop_id <- people$id[people$state == STATE_CONFIRMED &
187 people$die_moment == 0]
188 if ((npeop <- length(peop_id)) > 0) {
189 people$die_moment[peop_id] <- ifelse(
190 runif(npeop, 0, 1) < FATALITY_RATE, # 用均匀分布模拟确诊患者是否会死亡
191 people$confirmed_time + max(rnorm(npeop, DIE_TIME, DIE_SIGMA), 0), # 发病后确定死亡时刻
192 -1 # 逃过了死神的魔爪
193 )
194 }
195 # 如果患者已经确诊，且（世界时刻-确诊时刻）大于医院响应时间，
196 # 即医院准备好病床了，可以抬走了
197 peop_id <- people$id[people$state == STATE_CONFIRMED &
198 worldtime - people$confirmed_time >= HOSPITAL_RECEIVE_TIME]
199 if ((npeop <- length(peop_id)) > 0) {
200 if ((nbed_empty <- sum(hosp_beds$is_empty)) > 0) { # 有空余床位
201 nbed_use <- min(npeop, nbed_empty)
202 bed_id <- hosp_beds$id[hosp_beds$is_empty][1:nbed_use]
203 # 更新患者信息
204 peop_id2 <- sample(peop_id, nbed_use) # 这里是随机选择，理论上应该按症状轻重
205 people$x[peop_id2] <- hosp_beds$x[bed_id]
206 people$y[peop_id2] <- hosp_beds$y[bed_id]
207 people$state[peop_id2] <- STATE_FREEZE
208 people$freeze_time[peop_id2] <- worldtime
209 # 更新床位信息
210 hosp_beds$is_empty[bed_id] <- F
211 hosp_beds$person_id[bed_id] <- peop_id2
212 }
213 }
214 # TODO 需要确定一个变量用于治愈时长。
215 # 为了说明问题，暂时用一个正态分布模拟治愈时长并且假定治愈的人不会再被感染
216 peop_id <- people$id[people$state == STATE_FREEZE &
217 people$cured_moment == 0]
218 if ((npeop <- length(peop_id)) > 0) { # 正态分布模拟治愈时间
219 people$cured_moment[peop_id] <- people$freeze_time[peop_id] +
220 max(rnorm(npeop, CURED_TIME, CURED_SIGMA), 0)
221 }
222 peop_id <- people$id[people$state == STATE_FREEZE & people$cured_moment > 0 &
223 worldtime >= people$cured_moment]
224 if ((npeop <- length(peop_id)) > 0) { # 归还床位
225 people$state[peop_id] <- STATE_CURED
226 hosp_beds$is_empty[! hosp_beds$is_empty & hosp_beds$person_id %in% peop_id] <- T
227 people$x[peop_id] <- sapply(rnorm(npeop, CITY_CENTERX, PERSON_DIST_X_SIGMA),
228 format_coord, boundary = CITY_WIDTH) # (x, y) 为人群点坐标
229 people$y[peop_id] <- sapply(rnorm(npeop, CITY_CENTERY, PERSON_DIST_Y_SIGMA),
230 format_coord, boundary = CITY_HEIGHT)
231 people$tx[peop_id] <- rnorm(npeop, people$x[peop_id], PERSON_DIST_X_SIGMA)
232 people$ty[peop_id] <- rnorm(npeop, people$y[peop_id], PERSON_DIST_Y_SIGMA)
233 people$has_target[peop_id] <- T
234 people$is_arrived[peop_id] <- F
235 }
236 # 处理病死者
237 peop_id <- people$id[people$state %in% c(STATE_CONFIRMED, STATE_FREEZE) &
238 worldtime >= people$die_moment & people$die_moment > 0]
239 if (length(peop_id) > 0) { # 归还床位
240 people$state[peop_id] <- STATE_DEATH
241 hosp_beds$is_empty[! hosp_beds$is_empty & hosp_beds$person_id %in% peop_id] <- T
242 }
243 # 处理发病的潜伏期感染者
244 peop_id <- people$id[people$state == STATE_SHADOW &
245 worldtime >= people$confirmed_time]
246 if ((npeop <- length(peop_id)) > 0) {
247 people$state[peop_id] <- STATE_CONFIRMED # 潜伏者发病
248 }
249 # 处理未隔离者的移动问题
250 peop_id <- people$id[
251 ! people$state %in% c(STATE_FREEZE, STATE_DEATH) &
252 rnorm(CITY_PERSON_SIZE, MOVE_WISH_MU, MOVE_WISH_SIGMA) > 0] # 流动意愿
253 if ((npeop <- length(peop_id)) > 0) { # 正态分布模拟要移动到的目标点
254 pp_id <- peop_id[! people$has_target[peop_id] | people$is_arrived[peop_id]]
255 if ((npp <- length(pp_id)) > 0) {
256 people$tx[pp_id] <- rnorm(npp, people$tx[pp_id], PERSON_DIST_X_SIGMA)
257 people$ty[pp_id] <- rnorm(npp, people$ty[pp_id], PERSON_DIST_Y_SIGMA)
258 people$has_target[pp_id] <- T
259 people$is_arrived[pp_id] <- F
260 }
261 # 计算运动位移262 dx <- people$tx[peop_id] - people$x[peop_id]
263 dy <- people$ty[peop_id] - people$y[peop_id]
264 move_dist <- sqrt(dx ^ 2 + dy ^ 2)
265 people$is_arrived[peop_id][move_dist < 1] <- T # 判断是否到达目标点266 pp_id <- peop_id[move_dist >= 1]
267 if ((npp <- length(pp_id)) > 0) {
268 udx <- sign(dx[move_dist >= 1]) # x轴运动方向269 udy <- sign(dy[move_dist >= 1])
270 # 是否到了边界
271 pid_x <- (1:npp)[people$x[pp_id] + udx < 0 | people$x[pp_id] + udx > CITY_WIDTH]
272 pid_y <- (1:npp)[people$y[pp_id] + udy < 0 | people$y[pp_id] + udy > CITY_HEIGHT]
273 # 更新到了边界的点的信息
274 people$x[pp_id[pid_x]] <- people$x[pp_id[pid_x]] - udx[pid_x]
275 people$y[pp_id[pid_y]] <- people$y[pp_id[pid_y]] - udy[pid_y]
276 people$has_target[unique(c(pp_id[pid_x], pp_id[pid_y]))] <- F
277 # 更新没有到边界的点的信息278 people$x[pp_id[! pp_id %in% pid_x]] <- people$x[pp_id[! pp_id %in% pid_x]] +
279 udx[! pp_id %in% pid_x]
280 people$y[pp_id[! pp_id %in% pid_y]] <- people$y[pp_id[! pp_id %in% pid_y]] +
281 udy[! pp_id %in% pid_y]
282 }
283 }
284 # 处理健康人被感染的问题
285 # 通过一个随机幸运值和安全距离决定感染其他人286 normal_peop_id <- people$id[people$state == STATE_NORMAL]
287 other_peop_id <- people$id[! people$state %in% c(STATE_NORMAL, STATE_CURED)]
288 if (length(normal_peop_id) > 0) {
289 normal_other_dist <- apply(people[normal_peop_id, ], 1, get_min_dist,
290 peop = people[other_peop_id, ])
291 normal2other_id <- normal_peop_id[normal_other_dist < SAFE_DIST &
292 runif(length(normal_peop_id), 0, 1) < BROAD_RATE]
293 if ((n2other <- length(normal2other_id)) > 0) {
294 people$state[normal2other_id] <- STATE_SHADOW
295 people$infected_time[normal2other_id] <- worldtime
296 people$confirmed_time[normal2other_id] <- worldtime +
297 max(rnorm(n2other, SHADOW_TIME / 2, SHADOW_TIME_SIGMA), 0)
298 }
299 }
300
301 # 画出更新后的数据
302 npeople_confirmed <- sum(people$state >= STATE_CONFIRMED)
303 npeople_death <- sum(people$state == STATE_DEATH)
304 npeople_freeze <- sum(people$state == STATE_FREEZE)
305 npeople_shadow <- sum(people$state == STATE_SHADOW)
306 npeople_cured <- sum(people$state == STATE_CURED)
307 nbed_need <- npeople_confirmed - npeople_cured - npeople_death - BED_COUNT
308 nbed_need <- ifelse(nbed_need > 0, nbed_need, 0) # 不足病床数
309 # 疫情传播模拟散点图
310 dev.set(window_scatter)
311 plot(x = people$x, y = people$y, cex = .8, pch = 20, xlab = NA, ylab = NA,
312 xlim = c(5, max_plot_x), xaxt = "n", yaxt = "n", bty = "n", main = "疫情传播模拟",
313 sub = paste0("世界时间第 ", getday(worldtime), " 天"),
314 col = (people %>% left_join(person_color, by = "state") %>%
315 select(color))$color)
316 points(x = hosp_beds$x, y = hosp_beds$y, cex = .8, pch = 20,
317 col = (hosp_beds %>% left_join(bed_color, by = "is_empty") %>%
318 select(color))$color)
319 rect(HOSPITAL_X - BED_ROW_SPACE / 2, HOSPITAL_Y + 10 - BED_COLUMN_SPACE,
320 max(hosp_beds$x) + BED_ROW_SPACE / 2, max(hosp_beds$y + BED_COLUMN_SPACE))
321 legend(x = 150, y = -30, legend = person_color$label, col = person_color$color,
322 pch = 20, horiz = T, bty = "n", xpd = T)
323 # 人群变化模拟条形图
324 dev.set(window_hist)
325 bp_data <- c(npeople_death, npeople_cured, nbed_need, npeople_freeze,
326 npeople_confirmed, npeople_shadow)
327 bp <- barplot(bp_data, horiz = T, border = NA, col = bp_color,
328 xlim = c(0, CITY_PERSON_SIZE * 1), main = "人群变化模拟",
329 sub = paste0("世界时间第 ", getday(worldtime), " 天"))
330 abline(v = BED_COUNT, col = "gray", lty = 3)
331 abline(v = CITY_PERSON_SIZE, col = "gray", lty = 1)
332 text(x = -350, y = bp, labels = bp_labels, xpd = T)
333 text(x = bp_data + CITY_PERSON_SIZE / 15, y = bp, xpd = T,
334 labels = ifelse(bp_data > 0, bp_data, ""))
335 legend(x = 300, y = -.6, legend = c("总床位数", "城市总人口"), col = "gray",
336 lty = c(3, 1), bty = "n", horiz = T, xpd = T)
337
338 #更新世界时间
339 worldtime <- worldtime + 1
340 }相关文章