使用SIR模型对2019新型冠状病毒的疫情发展进行分析
使用SIR模型对2019新型冠状病毒的疫情发展进行分析
新型冠状病毒的确诊人数依旧在持续上升。在对传染病模型的研究上有很多模型比如：SI、SIS、SERS、SIR等，本文将利用SIR模型对这次新型冠状病毒的发展情况进行研究。
数据
数据本次数据比较简单可以看我之前文章爬取疫情数据，也可以直接直接手动输入。当然本次数据选取从一月份开始到2月12日，因为自从13日公布的确诊数据包含了临床数据，与之前的数据统计方式不一样因此没有加进去。那么先看下数据,在左边的图里，可以看到截止2月12日的确诊人数变化，右图是取完对数的变化并用线性模型拟合了一下，可以发现呈现出一种类似对数线性的关系。这表明该流行病处于指数阶段，尽管发病率略低于开始时。顺便说一句：一开始，很多人都没有惊慌。为什么？因为指数函数在开始时看起来是线性的。
相关python代码（python做个线性回归都要写个函数，所以接下来用R去建模）
import pandas as pd
import matplotlib.pyplot as plt
import math
import numpy as np
def linear_regression(x, y):
N = len(x)
sumx = sum(x)
sumy = sum(y)
sumx2 = sum(x ** 2)
sumxy = sum(x * y)
A = np.mat([[N, sumx], [sumx, sumx2]])
b = np.array([sumy, sumxy])
return np.linalg.solve(A, b)
data = pd.read_csv('data.csv',encoding='gb2312')
I = list(data['感染'])
N =1400000000
Day = []
logI = []
for i in range(len(I)):
Day.append(i+1)
logI.append(math.log(I[i]))
X1=np.array(Day)
Y1=np.array(logI)
a0, a1 = linear_regression(X1, Y1)
_Y1 = [a0 + a1 * x for x in Day]
ax1 = plt.subplot(1,2,1)
ax2 = plt.subplot(1,2,2)
plt.sca(ax1)
plt.scatter(Day,I, marker = 'x', s = 20,color='black')
plt.sca(ax2)
plt.scatter(Day,logI, marker = 'x', s = 20,color='black')
plt.yticks([])
plt.plot(Day, _Y1,color='red')SIR建模
SIR[1]模型比较简单，它将人群划分为三类人：健康但容易患病的人为易感人群(susceptible)，被感染的人(Infectious)和已康复的人(Recovered)，
为了模拟疫情爆发的动力学，我们需要三个微分方程:
这在R中很容易实现：
SIR <- function(time, state, parameters) {
par <- as.list(c(state, parameters))
with(par, { dS <- -beta/N * S * I
dI <- beta/N * S * I - gamma * I
dR <- gamma * I
list(c(dS, dI, dR))
})
}
接下来可以就交给R做参数估计和模型求解
library(deSolve)
Infected <- c(45, 62, 121, 198, 291, 440, 571, 830, 1287, 1975, 2744, 4515, 5974, 7711, 9692, 11791, 14380, 17205, 20440,24324,28018,31161,34546,37198,40171,42638,44653)
Day <- 0:(length(Infected)-1)
N <- 1400000000 #总人口
init <- c(S = N-Infected[1], I = Infected[1], R = 0)
RSS <- function(parameters) {
names(parameters) <- c("beta", "gamma")
out <- ode(y = init, times = Day, func = SIR, parms = parameters)
fit <- out[ , 3]
sum((Infected - fit)^2)
}
Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, 1)) #优化
Opt$message
Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par
## beta gamma
## 0.6746089 0.3253912
t <- 1:70 # 预测时间
fit <- data.frame(ode(y = init, times = t, func = SIR, parms = Opt_par))
col <- 1:3
matplot(fit$time, fit[ , 2:4], type = "l", xlab = "Day", ylab = "Number of subjects", lwd = 2, lty = 1, col = col)
matplot(fit$time, fit[ , 2:4], type = "l", xlab = "Day", ylab = "Number of subjects", lwd = 2, lty = 1, col = col, log = "y")
points(Day, Infected)
legend("bottomright", c("Susceptibles", "Infecteds", "Recovereds"), lty = 1, lwd = 2, col = col, inset = 0.05)
title("SIR model 2019-nCoV", outer = TRUE, line = -2)看下最后的结果，beta为0.6746089预测出来大概在两个月左右到达高峰，不过光凭简单的SIR模型估计不太好去准确预测，模型应该可以被进一步优化，同时从国家施行的各种管制措施，疫情应该得到了很好的控制。
也可以在R中计算
par(old)
R0 <- setNames(Opt_par["beta"] / Opt_par["gamma"], "R0")
R0
## R0
## 1.769682
fit[fit$I == max(fit$I), "I", drop = FALSE]
## I
## 61 157224962
最后
本次SIR建模分析的目的只是为了说明如何使用最简单的SIR模型，其结果依旧有很大的局限性。通过官方通报的部分病例来看，有些确诊病例的病毒潜伏期很长。尝试SEIR模型(Susceptible-Exposed-Infectious-Recovered)是否会更好可以进一步研究。最后，春天来了，也希望疫情能够尽快结束。
参考：
[1]SIR模型简介https://baike.baidu.com/item/SIR%E6%A8%A1%E5%9E%8B
[2]Novel coronavirus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions(https://doi.org/10.1101/2020.01.23.20018549)
腾讯云开发者
