使用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里面解微分方程可以用deSolve包，优化可以使用opt函数。方法采用的类似回归分析里的最小二乘，也就是使得预测与真实的之间的平方差最小：

接下来可以就交给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模型估计不太好去准确预测，模型应该可以被进一步优化，同时从国家施行的各种管制措施，疫情应该得到了很好的控制。

关于

也是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

可以看出为1.769682预计在60天左右感染人数达到峰值。注意到非典的约为0.85-3，埃博拉的为1.5-2.5。当然也不用害怕，本文建模过程中肯定有某些步骤可以被更好的优化，也一定有一些影响疫情发展的参数没有添加进来。不过在国外已经有研究者给出了较为准确的估计[2]，并且可以看出完全会影响到疫情的发展

最后

本次SIR建模分析的目的只是为了说明如何使用最简单的SIR模型，其结果依旧有很大的局限性。通过官方通报的部分病例来看，有些确诊病例的病毒潜伏期很长。尝试SEIR模型(Susceptible-Exposed-Infectious-Recovered)是否会更好可以进一步研究。最后，春天来了，也希望疫情能够尽快结束。

标签：SIR,plt,fit,模型,col,冠状病毒,2019,Day
来源： https://blog.51cto.com/15064626/2601427