一文秒懂病毒传播SIR模型 [附图表与代码]

marswriter

5 min readFeb 2, 2020

最近新型肺炎疫情引起各方关注，而且大家对数据似乎特别关注，都希望出现拐点，数字早日降下来。

(图片来自网络，各方支援！Be strong China 中国加油！🇨🇳🇨🇳🇨🇳)

有个著名的病毒传播学模型叫SIR模型，缩写自Susceptible(易感人群), Infected(感染人群), Recovered(免疫人群)。

敝人并不认为这个SIR模型符合本次新型肺炎病毒传播模型，因为模型假设所有人都会经历病毒并且能恢复，但新型肺炎显然是控制隔离来抑制，不然人类恐怕要灭种了。本文仅仅就数学模型本身做一些简单介绍。

模型很简单，根据前一天的情况推导出后一天的变化。

每日新发病数目跟易感人群(S)和感染人群(I)有关，设系数为a(即infect_rate)；
每日新治愈数目跟感染人群(I)有关，设系数为b(即recover_rate)；
总人数是一样的，就能得出每日感染人数(I)变化。

我们可以程序模拟一下，假设第一天，易感人群10万人，感染人群1千人，免疫人群仅20人，在不同的天感染率(infect_rate)和天治愈率(recover_rate)的情况下，看之后2个月的变化。

图表：Susceptible(易感人群), Infected(感染人群), Recovered(免疫人群)

低感染率，低治愈率，就是慢慢涨，2个月才过半

低感染率，高治愈率，似乎很快就被遏制了，希望疫情能这样早日被控制下去。

传染率略高，低治愈率，其实很快就所有人中招了，就是要慢慢恢复。

传染率略高，高治愈率，其实很快就所有人中招了，恢复很快。

高传染率，高治愈率，全军覆没，然后开始恢复。

附Python代码

import matplotlib.pyplot as plt
import pandas as pddef plot_line(input_df, xx, yy, label):
    plt.xlabel(xx)
    plt.plot(input_df[xx], input_df[yy], label=label)
    plt.legend()def sir_model(a, b, s0, i0, r0, days):
    """SIR model"""
    ss, ii, rr = [0]*days, [0]*days, [0]*days
    ss[0], ii[0], rr[0], total = s0, i0, r0, s0+i0+r0
    for t in range(1, days):
        delta1 = min(a*ss[t-1]*ii[t-1], ss[t-1])
        delta2 = min(ii[t-1], b*ii[t-1] - delta1)
        ss[t] = ss[t-1] - delta1  # susceptible decrease
        ii[t] = ii[t-1] - delta2
        rr[t] = total - ss[t] - ii[t]
    return ss, ii, rr
  
  
def sir_simulate(a, b):
    """simulate with a and b"""
    s_str, i_str, r_str = 'susceptible', 'infected', 'recovered'
    days = 60
    ss, ii, rr = sir_model(a=a, b=b, s0=100000, i0=1000, r0=20, days=days)trend = {
        'dt': [x for x in range(0, days)],
        s_str: ss,
        i_str: ii,
        r_str: rr,
    }df1 = pd.DataFrame(trend)
    for k in trend:
        if k == 'dt':
            continue
        plot_line(df1, 'dt', k, k)# main
arr_a = [0.000001, 0.00001, 0.0001]
arr_b = [0.01, 0.05, 0.1]
for a in arr_a:
    for b in arr_b:
        fig = plt.figure()
        title = 'infect_rate={}%, recover_rate={}%'.format(a*100, b*100)
        plt.title(title)
        sir_simulate(a, b)
        plt.show()

一文秒懂病毒传播SIR模型 [附图表与代码]

Written by marswriter

No responses yet