新冠病毒的SIR模型

化学课作业？

UPD：论文：
[pdf-embedder url=”https://www.yhzq-blog.cc/wp-content/uploads/2020/03/Covid-19.pdf”]

化学课的期末作业是对新冠病毒进行数学建模。。也不需要特别精确，就用最简单的SIR模型吧。

S代表易感者，I表示感染者，R表示恢复者+死亡者。易感者会变成感染者，感染者会变成恢复者或者死亡者（这两种没有区别，都不会再传染或者被传染。那么这可以用一种化学反应表示：

$$S + I \mathop\rightarrow^r 2I \\ I \mathop\rightarrow^a R$$

第一个是感染的过程，反应速率为r，I相当于一种自催化剂。第二个是恢复和死亡的过程，速率为a。根据rate law，我们能列出如下ODE

$${dS\over dt} = -rSI \\ {dI\over dt} = rSI – aI \\ {dR\over dt} = aI$$

这就是SIR模型了。这个模型的一些假设：

1.人口在一段时间内保持不变。
2.得了传染病的人会有一定的传染性。
3.。。。其实还有很多，想想就不合理。比如潜伏期啥的

在这里面还有一个重要的参数叫基本再生数$R_0$。它表示在发病初期，一个病人在平均患病期能传染的人数。所以$R_0$和1的大小关系取决了这个病是否能传播。通过这个ODE：
$${dR\over dt} = aI$$
可以解得
$$R(t) ={e^{aIt} + C}\\ Since\ R(0) = 0, R(t) = e^{aIt}$$
我们让$I = 1$，且不考虑I的变化，那么这个函数就代表了这个一个患者在t转换成死亡或者痊愈的概率。同理可得，一个患者在t仍然存活的概率为$e^{-at}$。那么，在t时刻新增的患者（由这一个患者传染的）为$rS_0e^{-at}$，所以
$$R_0=\int_0^\infty rS_0e^{-at}dt={rS_0\over a}$$

首先在中国健康部爬下数据。我决定使用湖北的数据，因为假设2，最好使用一个小区域的数据。而且湖北是病毒爆发的地方。湖北的数据：

现在我们有三个变量不知道：$S_0,r,a$, 分别是易感人群的初始数量和两个速率。首先，我们知道$1/a$应该为恢复的天数。查了些资料发现大概是25天，就先假设$a=0.04$。剩下两个利用双重循环暴力试下。利用R^2做最优。

解ODE的函数：

def get_a(n, S, r, a):
    I = 494
    R = 0
    dt = 0.001
    dS = 0
    dI = 0
    dR = 0
    aS = [S]
    aI = [I]
    aR = [R]
    cont = 0
    for i in range(int(n / dt)):
        S += dS
        I += dI
        R += dR
        if (cont == 1 / dt):
            aS.append(S)
            aI.append(I)
            aR.append(R)
            cont = 0
        cont += 1
        dS = -1 * r * S * I * dt
        dI = (r * S * I - a * I) * dt
        dR = a * I * dt
    return np.asarray(aS), np.asarray(aI), np.asarray(aR)

下面是找最优的

def get_error(aI, aR, aS, S):
    mI = sum((aI - tI) ** 2) / len(tI)
    mR = sum((aR - tR) ** 2) / len(tR)
    mS = sum((aS - (S - tI - tR)) ** 2) / len(aS)
    rI = mI / np.var(tI)
    rR = mR / np.var(tR)
    rS = mS / np.var(S - tI - tR)
    return (rI + rR + rS) / 3
def get_max(n, a):
    bS, bI, bR = get_a(n, 60000, 1e-5, a)
    bE = get_error(bI, bR, bS, 60000)
    ers = []
    for S in range(68000, 100000, 1000):
        cS, cI, cR = get_a(n, S, 1e-5, a)
        cE = get_error(cI, cR, cS, S)
        po = -1
        for r in np.arange(1e-7, 1e-5, 5e-7):
            aS, aI, aR = get_a(n, S, r, a)
            tE = get_error(aI, aR, aS, S)
            if tE < cE:
                cS, cI, cR = aS, aI, aR
                cE = tE
                po = r
        ers.append(cE)
        if cE < bE:
            bS, bI, bR = cS, cI, cR
            bE = cE
            print(po, S, 1 - bE, po * S / a)
    return bS, bI, bR, ers

返回的ers可以看出现在的人数区间是否有最好的。比如画出来是这：

能看出来已经经过极值点了，这个区间就是可行的。找出比较小的r和S的区间后，再加入a的变化。

mins = 1e15
for a in np.arange(0.035, 0.042, 0.001):
    bS, bI, bR, cE = get_max(54, a)
    if np.min(cE) < mins:
        mins = np.min(cE)
        aS, aI, aR = bS, bI, bR
        print(a, 1 - mins)

然后就得出了这样的结果：