前言
在文章开始前，我得引用一下《关于新冠肺炎的一切｜回形针》的一条数据，这个数据的最终来源应该是《公共卫生新研究：中国流感相关呼吸死亡负担严重，特别是老年人群 | 柳叶刀》：
过去几年，中国平均每年有 8.8 万人死于流感引发的呼吸系统疾病
中国的舆论也是蛮有趣的，2020 年 2 月 22 日很多中国网民都在微博和今日头条底下说新冠病毒是先从美国爆发（因为美国正在爆发流感，还死了几万人，中国网民认为美国的流感就是新冠病毒引起的），然后通过军运会传到武汉，部分微博舆论截图如下所示：
部分今日头条舆论截图如下所示：
如果您仔细想想，这种舆论本身就是矛盾的：
- 如果是从美国开始爆发，那么最先遭殃的应该是北上广深这些大城市，因为这些城市国际航班比较多，而不是先从武汉开始。
- 如果是美国士兵通过军运会带来武汉的，那么新冠病毒不可能会在美国本土爆发吧？自研的生化武器先在美国本土做测试再带来武汉？杀敌一千，自损八百？而且相比于大城市的国际航班，一个国家来不了多少军人的。
病毒性感冒和病毒性肺炎还是不一样的，病毒性感冒大部分人都可以自愈，而病毒性肺炎是肺部感染，这往往容易造成呼吸困难（类似于溺水），从而引发其他并发症，这就需要占用大量的医疗资源（至少呼吸困难了要吸氧）。
当然了，也不是没有明理的人，部分微博舆论截图如下所示：
可能有的人会说 “ 中国平均每年有 8.8 万人死于流感引发的呼吸系统疾病 ” 是估计的，不是层层上报的准确的数据。这个数据确实是估计的，《公共卫生新研究：中国流感相关呼吸死亡负担严重，特别是老年人群 | 柳叶刀》里面有提到：
流感的临床特征不特异，实验室检测是确诊流感的唯一手段，但目前仅采集少数上呼吸道感染症状的病人标本进行实验室确诊。因此，基于实验室确诊的流感死亡病例仅是其死亡负担的冰山一角，越来越多的研究采用基于监测数据的统计建模方法来估计流感的死亡负担。
但同时也提到 H1N1 猪流感病毒也是有在中国传播开来的，这篇文章估计的是 “ 甲型H1N1流感大流行后全国和各省流感的死亡负担 ” ：
前期已有研究阐明了2009年甲型H1N1流感大流行前中国部分地区的流感死亡负担[4]，以及甲型H1N1流感大流行前和大流行期间全国流感的死亡负担[5]，但目前尚无研究报道甲型H1N1流感大流行后全国和各省流感的死亡负担。针对上述科学问题，本研究通过统计建模的方法，精确估计了2010-2011至2014-2015年间，全国和各省流感相关的超额呼吸死亡率，为中国制定国家级及区域的流感防控措施提供重要科学依据。
而且美国疾病控制与预防中心 CDC 的数据也是预估的：
2019到2020流感季（从2019年10月1日到2020年2月15日），预计至少2900万人生病，至少28万人住院，至少16000人死于流感和肺炎等并发症。
同时，钟南山院士在 2020 年 2 月 23 日也提到，中国也有美国类似的情况，即中国也有流感，把两者（流感和新冠肺炎）鉴别开来是当务之急，如果流感和新冠肺炎都是一回事儿那还区分个啥？该新闻来自《钟南山：当务之急要鉴别流感和新冠肺炎 | 广州日报》。
2020 年 2月 29 日一段来自台湾的节目《这！不是新闻》的视频片段在微博火了，这个节目是 EBC 东森财经新闻播出的时事新闻、财经谈话类节目。以下是部分反驳该台湾节目的微博博主（ ZHANKAI9426 、果壳、骆明微博和 fengfeixue0219 ）：
可能有的人会说，上面这些不是官方媒体说的，那我们来看看官方媒体《新冠病毒源头在美国？专家解读来了 | 经济日报》，其中有一段：
据此，并结合病患发病时间记录和种群扩张时间，论文研究团队推断：华南海鲜市场的新型冠状病毒是从其它地方传入进来，在市场中发生快速传播蔓延到市场之外。换句话说，该论文认为，华南海鲜市场不是病毒发源地。
但同时也提到：
因此，基于这篇论文我们无法得出新冠病毒源头在美国的结论。
“根据基因组分析，新冠病毒的源头并不在华南海鲜市场。其源头究竟在哪，还需要相关研究者的继续追踪，但根据目前的相关研究，应该还在武汉。”赵序茅称。
华南海鲜市场不是病毒发源地这个观点 Ricky 还是认可的，但最终来源还需要做大量的追踪调查，目前很难肯定地说新冠病毒的源头就是美国，或者新冠病毒的源头就不是美国。
2003 年 SARS 在中国广东爆发，石正丽相关团队在经过 13 年的追踪最终发现 SARS 病毒起源于云南某山洞里的蝙蝠，该消息来自《这些野生动物的病毒怎么就到了人类社会？| 石正丽 一席第600位讲者》：
一直到了2011年，我们有了新的线索：我们在云南一个洞里面分离并检测到了和SARS病毒高度同源的，可以说是SARS的直系亲属的一株蝙蝠SARS样冠状病毒。除了进化关系近，在功能上面也非常相近。
……
经过13年的追踪，我们最终确定了SARS病毒的蝙蝠起源。但是当年SARS爆发于广东，我们发现的病毒是在云南，它是怎么过去的呢？
其实SARS暴发前，我国很多地方养殖的果子狸，绝大部分销售到了广东。我们推测，蝙蝠SARS样冠状病毒在偶然的情况下感染了云南养殖场的果子狸，感染了病毒的果子狸随后又被贩卖到了广东。病毒进一步在市场上的果子狸中传播，不断变异，最终产生一个传播性极强的SARS病毒，感染了人类。
既然蝙蝠才是SARS的元凶，果子狸是不是被冤枉了？也不是。虽然SARS的根源不在果子狸，但毕竟它是人感染SARS病毒的直接来源。
Ricky 想说的是，这些阴谋论患者还是省省吧，有这功夫不如来学习一下传染病模型。毕竟恐惧往往来源于未知，消除恐惧的最好办法就是面对恐惧。看完本篇文章以后您或许会对流感或新冠肺炎的传染过程有了更加深刻的理解。【没错！强行安利！】
正文开始
2020 年 1 月 20 日前生活依旧，直到钟南山院士公开说肯定是人传人以及 1 月 23 日武汉宣布封城的时候，Ricky 才预感到事态的严重。
一月底，当确诊人数还在三位数至四位数出头的时候，Ricky 在微信里看到一些非官方的小视频，有护士说至少 1 万人感染，Ricky 当时都有点不相信，毕竟 2003 年 SARS 官方统计的感染人数都没有那么多（截至 2003 年 8 月 16 日患病人数 8422 例，死亡人数 919 例，数据来源：百度百科）。当时正好有看到一些传染病模型的介绍，再加上有人做了些预测（如：《简单算算，你宅在家里究竟能为抗击肺炎疫情做出多大贡献？| bilibili 》），所以 Ricky 打算也自己试着做做看。
如有不足之处，欢迎在下方留言并批评指正。
模型一 ：SI 模型
我们先从最简单的 SI 模型开始说起，SI 模型把人群分为两种：
- S ：Susceptibles ，易感者（也就是健康者或健康人群）
- I ：The Infected ，感染者
如果我们假设一个区域内的总人数是 ，那么 。
如果我们再假设每天有 个感染者到处溜达，每个感染者每天接触 个人，有 的概率会将病毒传染给易感者，那么其微分方程如下所示：
看到这里我相信很多人都看懵圈了，别急，我一个个解释：
- 所有感染者每天接触到的总人数是 （注意：接触到的人中既有感染者又有易感者）；
- 因为所有易感者在区域内的占比是 ，所以所有感染者每天接触到的易感者是 ；
- 那么每天新增的感染者就是 ，可得 ；
- 相对的，每天减少的易感者就是 ，可得 。
我们举个例子：如果我们假设一个区域内的总人数是 1000 人，我们再假设今天有 5 个感染者到处溜达，每个感染者每天接触 50 个人，有 1 % 的概率会将病毒传染给易感者。那么：
- 今天所有感染者接触到的总人数是 5 x 50 = 250 ，这 250 个人中什么人都有，因为也有可能感染者还会接触到另一个感染者；
- 因为所有易感者在区域内的占比是 ( 1000 – 5 ) / 1000 = 995 / 1000 ，这其实就是感染者遇到易感者的概率，所以今天所有感染者接触到的易感者是 5 x 50 x 995 / 1000 = 248.75 人；
- 因为感染者有 1 % 的概率会将病毒传染给易感者，所以这 1000 个人中今天新增的感染者就是 5 x 50 x 1 % x 995 / 1000 = 2.4875 人；
- 995 个易感者中今天有 2.4875 人感染，那么易感者人群就得减掉这 2.4875 人，即今天的易感者还剩下 992.5125 人。
有了这些公式，我们就可以进行仿真了，上述式子实际上就是一个马尔可夫链。以下解释转自《马尔可夫与隐马尔可夫模型 | 知乎》：
马尔可夫链（英语：Markov chain），又称离散时间马尔可夫链（discrete-time Markov chain，缩写为DTMC），因俄国数学家安德烈·马尔可夫（俄语：Андрей Андреевич Марков）得名，为状态空间中经过从一个状态到另一个状态的转换的随机过程。该过程要求具备“无记忆”的性质：下一状态的概率分布只能由当前状态决定，在时间序列中它前面的事件均与之无关。这种特定类型的“无记忆性”称作马尔可夫性质。马尔科夫链作为实际过程的统计模型具有许多应用。
所以，第 天的人数我们可以表示为：
每天新增的感染者：
使用 Python 3.8 写的示例代码如下所示：
def SI_model(N = 10000, I = 1, r = 10, b = 0.01, days = 200): ''' I 感染者人数（ The Infected ） S 易感者人数（ Susceptibles ） N 区域内总人口，N = I + S S / N 【易感者人数】占【区域内总人口】的比例 r 每个感染者每天接触的人数 b 感染者将病毒传染给易感者的概率 days 是天数 ''' S = N - I I_list = [I] S_list = [S] for i in range(1, days): I_perday = r * b * I_list[i - 1] * S_list[i - 1] / N I_list.append(I_list[i - 1] + I_perday) S_list.append(S_list[i - 1] - I_perday)
我们建立一个 SI 模型示例看看，设置的参数如下所示：
- 区域内的总人口
- 感染者人数
- 每个感染者每天接触的人数
- 感染者将病毒传染给易感者的概率
- 模型需要模拟的总天数为
我们使用 Python 3.8 和 ReportLab 进行模拟：
上图中：The Infected 是感染者人数，Susceptibles 是易感者人数，Total Population 是总人口。
《生化危机》里变丧尸就是这个模型，感染者不可恢复，现实生活中就是艾滋病。
模型二：SIS 模型
我们在 SI 模型的基础上增加复杂度，有的人治好了但还是会反复感染，类似流感。SIS 模型在 SI 模型的基础上多了一个参数：
- ：每天感染者被治愈的人数占感染者总数的比例（日治愈率）
也就是说易感者 会变成感染者 ，感染者 又会变回易感者 ，这就是这个模型名称的由来（ SIS ）。
该 SIS 模型的微分方程如下所示：
这个模型应该很容易理解：
- 因为每天总的感染者人数是 ，所以每天被治愈的人数是 ；
- 因为每天新增的感染者是 ，每天被治愈的人数是 ，所以每天感染者的变化人数是 （要从感染者中减去治愈的人数）；
- 同理，每天易感者的变化人数就是 。
所以，第 天的人数我们可以表示为：
每天新增的感染者：
每天被治愈的人数：
使用 Python 3.8 写的示例代码如下所示：
def SIS_model(N = 10000, I = 1, r = 10, b = 0.01, y = 0.02, days = 200): ''' I 感染者人数（ The Infected ） S 易感者人数（ Susceptibles ） N 区域内总人口，N = I + S S / N 【易感者人数】占【区域内总人口】的比例 r 每个感染者每天接触的人数 b 感染者将病毒传染给易感者的概率 days 是天数 y 每天感染者被治愈的人数占感染者总数的比例 ''' S = N - I I_list = [I] S_list = [S] for i in range(1, days): I_perday = r * b * I_list[i - 1] * S_list[i - 1] / N S_perday = y * I_list[i - 1] I_list.append(I_list[i - 1] + I_perday - S_perday) S_list.append(S_list[i - 1] - I_perday + S_perday)
我们建立一个 SIS 模型示例看看，设置的参数如下所示：
- 区域内的总人口
- 感染者人数
- 每个感染者每天接触的人数
- 感染者将病毒传染给易感者的概率
- 每天感染者被治愈的人数占感染者总数的比例
- 模型需要模拟的总天数为
我们使用 Python 3.8 和 ReportLab 进行模拟：
上图中：The Infected 是感染者人数，Susceptibles 是易感者人数，Total Population 是总人口。
这个模型适用于日常生活中的流感，我们可以看到感染者和易感者最终形成了一个稳定的动态平衡。
模型三：SIR 模型
我们继续在 SIS 模型的基础上增加复杂度，我们假设人在康复以后产生了抗体就不会再得病。此时，我们需要在 SIR 模型中引入一个新的人群：
- R ：The Recovered ，康复者
也就是说易感者 会变成感染者 ，感染者 又会变成康复者 ，这就是这个模型名称的由来（ SIR ）。
如果我们假设一个区域内的总人数是 ，那么 。
该 SIR 模型的微分方程如下所示：
这个模型应该也很容易理解：
- 由 SIS 模型可知，每天感染者的变化人数是 ，这个没有变化；
- 在 SIS 模型中，每天被治愈的人数 是需要重新转化为易感者的，所以每天易感者的变化人数是 ；
- 而在 SIR 模型中，每天被治愈的人数 不会转化为易感者，而是转化为康复者，所以每天康复者的变化人数是 ，每天易感者的变化人数是 。
所以，第 天的人数我们可以表示为：
每天新增的感染者：
每天新增的康复者：
使用 Python 3.8 写的示例代码如下所示：
def SIR_model(N = 10000, I = 1, r = 10, b = 0.05, y = 0.1, days = 100): ''' I 感染者人数（ The Infected ） S 易感者人数（ Susceptibles ） R 康复者人数（但人在康复以后产生了抗体就不会再得病，The Recovered ） N 区域内总人口，N = I + S + R S / N 【易感者人数】占【区域内总人口】的比例 r 每个感染者每天接触的人数 b 感染者将病毒传染给易感者的概率 days 是天数 y 每天感染者被治愈的人数占感染者总数的比例 ''' S = N - I I_list = [I] S_list = [S] R_list = [0] for i in range(1, days): I_perday = r * b * I_list[i - 1] * S_list[i - 1] / N R_perday = y * I_list[i - 1] I_list.append(I_list[i - 1] + I_perday - R_perday) S_list.append(S_list[i - 1] - I_perday) R_list.append(R_list[i - 1] + R_perday)
我们建立一个 SIR 模型示例看看，设置的参数如下所示：
- 区域内的总人口
- 感染者人数
- 每个感染者每天接触的人数
- 感染者将病毒传染给易感者的概率
- 每天感染者被治愈的人数占感染者总数的比例
- 模型需要模拟的总天数为
我们使用 Python 3.8 和 ReportLab 进行模拟：
上图中：The Infected 是感染者人数，Susceptibles 是易感者人数，The Recovered 是康复者人数，Total Population 是总人口。
适合这个模型的疾病有：麻疹和腮腺炎。由上图可知，在该疫情下第 30 天左右感染人数会达到一个峰值，也就是大约 5000 人；但随着康复人数的增加，拥有抗体的人数也随即增加，病情就将得到控制。
模型四：SEIR 模型
我们继续在 SIR 模型的基础上增加复杂度，我们假设易感者在一开始会经历一段潜伏期，一段时间之后才会出现症状。此时，我们需要在 SIER 模型中引入一个新的人群：
- E ：The Exposed ，潜伏者
也就是说易感者 会先变成潜伏者 ，然后再变成感染者 ，感染者 又会变成康复者 ，这就是这个模型名称的由来（ SEIR ）。
如果我们假设一个区域内的总人数是 ，那么 。
同时 SEIR 模型比 SIR 模型多了一个参数：
- ：潜伏者每天按照概率 转化为感染者（ 在数学意义上等同于潜伏期时长的倒数）
该 SEIR 模型的微分方程如下所示：
这个模型的理解也不是很难：
- 因为每天总的潜伏者人数是 ，所以每天从潜伏者转化为感染者的人数是 ；
- 这里需要注意的是 是感染者把病毒传染给易感者的概率，但在把病毒传染给易感者后，易感者应该转化为潜伏者而不是感染者，所以 是每天易感者转化为潜伏者的人数；
- 综上，每天易感者的变化人数是 ，这部分减去的易感者是每天从易感者转化为潜伏者的人数（转出）；
- 每天潜伏者的变化人数是 ，除了每天会有易感者转化为潜伏者以外（转入），还会有潜伏者转化为感染者（转出，所以要 ）；
- 每天感染者的变化人数是 ，除了每天会有潜伏者转化为感染者以外（转入），还会有感染者转化为康复者（转出，所以要 ）；
- 每天康复者的变化人数是 ，康复者自然是从感染者那边转化过来的。
所以，第 天的人数我们可以表示为：
每天从易感者转化为潜伏者的人数：
每天从潜伏者转化为感染者的人数：
每天新增的康复者：
使用 Python 3.8 写的示例代码如下所示：
def SEIR_model(N = 10000, E = 0, I = 1, r = 20, b = 0.03, y = 0.1, a = 0.1, days = 140): ''' I 感染者人数（ The Infected ） S 易感者人数（ Susceptibles ） E 潜伏者人数（ 易感者一开始会经历一段潜伏期，一段时间之后才会出现症状，The Exposed ） a 潜伏者按照概率 a 转化为感染者（ a 在数学意义上等同于潜伏期时长的倒数） R 康复者人数（但人在康复以后产生了抗体就不会再得病，The Recovered ） N 区域内总人口，N = I + S + R + E S / N 【易感者人数】占【区域内总人口】的比例 r 每个感染者每天接触的人数 b 感染者将病毒传染给易感者的概率 days 是天数 y 每天感染者被治愈的人数占感染者总数的比例 ''' S = N - I - E I_list = [I] S_list = [S] E_list = [E] R_list = [0] for i in range(1, days): E_perday = r * b * I_list[i - 1] * S_list[i - 1] / N I_perday = a * E_list[i - 1] R_perday = y * I_list[i - 1] I_list.append(I_list[i - 1] + I_perday - R_perday) S_list.append(S_list[i - 1] - E_perday) R_list.append(R_list[i - 1] + R_perday) E_list.append(E_list[i - 1] + E_perday - I_perday)
我们建立一个 SEIR 模型示例看看，设置的参数如下所示：
- 区域内的总人口
- 感染者人数
- 潜伏者人数
- 每个感染者每天接触的人数
- 感染者将病毒传染给易感者的概率
- 每天感染者被治愈的人数占感染者总数的比例
- 潜伏者每天转化为感染者的概率
- 模型需要模拟的总天数为
我们使用 Python 3.8 和 ReportLab 进行模拟：
上图中：The Infected 是感染者人数，Susceptibles 是易感者人数，The Recovered 是康复者人数，The Exposed 是潜伏者人数，Total Population 是总人口。
这个模型适用于 SARS 病毒。但这还不是新冠病毒所使用的模型，因为新冠肺炎的潜伏期也能感染人，所以还需要另一个 SEIR 模型。
模型五：SEIR 2.0 模型
现在已经非常接近新冠病毒了，但我们知道这次新冠肺炎在潜伏期也具有传染性（即不光是感染者能将易感者转化为潜伏者，潜伏者也能将易感者转化为潜伏者），因此我们需要在 SEIR 模型的基础上多引入两个参数：
- ：每个潜伏者每天接触的人数
- ：潜伏者把病毒传染给易感者的概率
该 SEIR 2.0 模型的微分方程如下所示：
这个模型跟 SEIR 模型基本是一致的：
- 所有潜伏者每天接触到的总人数是 （注意：接触到的人中潜伏者、感染者和易感者都有）；
- 因为所有易感者在区域内的占比是 ，所以所有潜伏者每天接触到的易感者是 ；
- 那么每天由潜伏者把易感者转化为潜伏者的人数就是 。
所以，第 天的人数我们可以表示为：
每天由感染者把易感者转化为潜伏者的人数：
每天由潜伏者把易感者转化为潜伏者的人数：
每天从潜伏者转化为感染者的人数：
每天新增的康复者：
使用 Python 3.8 写的示例代码如下所示：
def SEIR_2_model(N = 10000, E = 0, I = 1, r = 20, r2 = 20, b = 0.03, b2 = 0.03, y = 0.1, a = 0.1, days = 140): ''' I 感染者人数（ The Infected ） S 易感者人数（ Susceptibles ） E 潜伏者人数（ 易感者一开始会经历一段潜伏期，一段时间之后才会出现症状，The Exposed ） a 潜伏者按照概率 a 转化为感染者（ a 在数学意义上等同于潜伏期时长的倒数） R 康复者人数（但人在康复以后产生了抗体就不会再得病，The Recovered ） N 区域内总人口，N = I + S + R + E S / N 【易感者人数】占【区域内总人口】的比例 r 每个感染者每天接触的人数 b 感染者将病毒传染给易感者的概率 r2 每个潜伏者每天接触的人数 b2 潜伏者将病毒传染给易感者的概率 days 是天数 y 每天感染者被治愈的人数占感染者总数的比例 ''' S = N - I - E I_list = [I] S_list = [S] E_list = [E] R_list = [0] for i in range(1, days): E_perday = r * b * I_list[i - 1] * S_list[i - 1] / N E_perday_for_E = r2 * b2 * E_list[i - 1] * S_list[i - 1] / N I_perday = a * E_list[i - 1] R_perday = y * I_list[i - 1] I_list.append(I_list[i - 1] + I_perday - R_perday) S_list.append(S_list[i - 1] - E_perday - E_perday_for_E) R_list.append(R_list[i - 1] + R_perday) E_list.append(E_list[i - 1] + E_perday + E_perday_for_E - I_perday)
我们建立一个 SEIR 2.0 模型示例看看，设置的参数如下所示：
- 区域内的总人口
- 感染者人数
- 潜伏者人数
- 每个感染者每天接触的人数
- 每个潜伏者每天接触的人数
- 感染者将病毒传染给易感者的概率
- 潜伏者将病毒传染给易感者的概率
- 每天感染者被治愈的人数占感染者总数的比例
- 潜伏者每天转化为感染者的概率
- 模型需要模拟的总天数为
我们使用 Python 3.8 和 ReportLab 进行模拟：
上图中：The Infected 是感染者人数，Susceptibles 是易感者人数，The Recovered 是康复者人数，The Exposed 是潜伏者人数，Total Population 是总人口。
这个模型基本上适用于新冠病毒了。我们能发现，因为潜伏期也具备传染性，所以此次疫情比 SEIR 模型示例中描述的疫情更严重了。
在这里非常感谢各位看官，因为这么冗长的文章您们还能耐着性子看到这里。那么文章是否要结束了呢？别急还有两个模型和对新冠肺炎的简单预测 ……
模型六：SEIR 3.0 模型
SEIR 2.0 模型都说到新冠肺炎了，SEIR 3.0 模型又是什么鬼？因为 SEIR 2.0 模型我们没有考虑到人为的因素，也就是封城（甚至是戒严，可参考《中华人民共和国戒严法 | 百度百科》）。相对于其他一些国家来说可能就是宣布国家进入紧急状态，从而进行 “ 战时 ” 管理或者战时管理。
SEIR 3.0 模型同样也是基于上一个模型 —— SEIR 2.0 模型发展而来的，SEIR 3.0 模型在 SEIR 2.0 模型的基础上新增了三个参数：
- ：是从第几天开始戒严
- ：是戒严后的
- ：是戒严后的
微分方程这些与 SEIR 2.0 模型一致，主要的不同在代码逻辑上，使用 Python 3.8 写的示例代码如下所示：
def SEIR_3_model(N = 10000, E = 0, I = 1, r = 20, r2 = 20, re = 5, r2e = 5, martial_law_day = 10, b = 0.03, b2 = 0.03, y = 0.1, a = 0.1, days = 140):
'''
I 感染者人数（ The Infected ）
S 易感者人数（ Susceptibles ）
E 潜伏者人数（ 易感者一开始会经历一段潜伏期，一段时间之后才会出现症状，The Exposed ）
a 潜伏者按照概率 a 转化为感染者（ a 在数学意义上等同于潜伏期时长的倒数）
R 康复者人数（但人在康复以后产生了抗体就不会再得病，The Recovered ）
N 区域内总人口，N = I + S + R + E
S / N 【易感者人数】占【区域内总人口】的比例
r 每个感染者每天接触的人数
re 是戒严后的 r
b 感染者将病毒传染给易感者的概率
r2 每个潜伏者每天接触的人数
r2e 是戒严后的 r2
b2 潜伏者将病毒传染给易感者的概率
days 是天数
martial_law_day 是从第几天开始戒严
y 每天感染者被治愈的人数占感染者总数的比例
'''
S = N - I - E
I_list = [I]
S_list = [S]
E_list = [E]
R_list = [0]
for i in range(1, days):
if (i + 1) == martial_law_day:
r, r2 = re, r2e
E_perday = r * b * I_list[i - 1] * S_list[i - 1] / N
E_perday_for_E = r2 * b2 * E_list[i - 1] * S_list[i - 1] / N
I_perday = a * E_list[i - 1]
R_perday = y * I_list[i - 1]
I_list.append(I_list[i - 1] + I_perday - R_perday)
S_list.append(S_list[i - 1] - E_perday - E_perday_for_E)
R_list.append(R_list[i - 1] + R_perday)
E_list.append(E_list[i - 1] + E_perday + E_perday_for_E - I_perday)
我们建立一个 SEIR 3.0 模型示例看看，设置的参数如下所示：
- 区域内的总人口
- 感染者人数
- 潜伏者人数
- 每个感染者每天接触的人数
- 是戒严后的 ，
- 每个潜伏者每天接触的人数
- 是戒严后的 ，
- 从第几天开始戒严
- 感染者将病毒传染给易感者的概率
- 潜伏者将病毒传染给易感者的概率
- 每天感染者被治愈的人数占感染者总数的比例
- 潜伏者每天转化为感染者的概率
- 模型需要模拟的总天数为
我们使用 Python 3.8 和 ReportLab 进行模拟：
上图中：The Infected 是感染者人数，Susceptibles 是易感者人数，The Recovered 是康复者人数，The Exposed 是潜伏者人数，Total Population 是总人口。
我们可以发现，相比于 SEIR 2.0 模型示例，在 SEIR 3.0 模型示例中疫情得到了极大的控制，这包含两个方面：
- 疾病高发期的时间向后延长
- 潜伏者和感染者高峰人数的大幅下降
这将为抗击新冠病毒争取到更多的时间，同时感染者高峰人数的下降也将大幅降低对医疗资源乃至社会资源的消耗。简而言之，封城（或者相对其他一些国家来说是宣布国家进入紧急状态）确实是行之有效的办法，而且效果明显。
再坚持一下，Ricky 也是快写傻了，还剩最后一个模型。
模型七：SEIR(D) 4.0 模型
那 SEIR(D) 4.0 模型又是什么鬼？这次我们打算考虑进来一个外部因素 —— 戴口罩，同时尝试引入一个新的人群：
- D ：Death ，病死者
加入 D 的目的是为了预测累计死亡人数，部分感染者 会转化为病死者 。
如果我们假设一个区域内的总人数是 ，那么 。
SEIR(D) 4.0 模型同样也是基于上一个模型 —— SEIR 3.0 模型发展而来的，SEIR(D) 4.0 模型在 SEIR 3.0 模型的基础上新增了两个参数：
- ：相当于是口罩的防护率，戒严开始后： ， （医用口罩的防护率一般是 ）
- ：每天感染者死亡的人数占感染者总数的比例（日死亡率）
该 SEIR(D) 4.0 模型的微分方程如下所示：
如果前面几个模型您都能看懂，那么这个模型似乎就不需要解释了：
- 每天感染者死亡的人数是 ，这部分人会从每天感染者的变化人数 中减去，然后加入到 中。
所以，第 天的人数我们可以表示为：
每天由感染者把易感者转化为潜伏者的人数：
每天由潜伏者把易感者转化为潜伏者的人数：
每天从潜伏者转化为感染者的人数：
每天新增的康复者：
每天新增的病死者：
使用 Python 3.8 写的示例代码如下所示：
def __SEIRD_4_model(N, E, I, r, r2, re, r2e, p, d, martial_law_day, b, b2, y, a, days): ''' I 感染者人数（ The Infected ） S 易感者人数（ Susceptibles ） E 潜伏者人数（ 易感者一开始会经历一段潜伏期，一段时间之后才会出现症状，The Exposed ） a 潜伏者按照概率 a 转化为感染者（ a 在数学意义上等同于潜伏期时长的倒数） R 康复者人数（但人在康复以后产生了抗体就不会再得病，The Recovered ） N 区域内总人口，N = I + S + R + E S / N 【易感者人数】占【区域内总人口】的比例 r 每个感染者每天接触的人数 re 是戒严后的 r b 感染者将病毒传染给易感者的概率 r2 每个潜伏者每天接触的人数 r2e 是戒严后的 r2 b2 潜伏者将病毒传染给易感者的概率 p 相当于是口罩的防护率，戒严开始后，b = b * p ， b2 = b2 * p days 是天数 martial_law_day 是从第几天开始戒严 y 每天感染者被治愈的人数占感染者总数的比例 d 每天感染者死亡的人数占感染者总数的比例 ''' S = N - I - E I_list = [I] S_list = [S] E_list = [E] R_list = [0] D_list = [0] C_list = [0] # Confirmed list ，为了记录累计感染者人数 for i in range(1, days): if (i + 1) == martial_law_day: r, r2, b, b2 = re, r2e, b * p, b2 * p E_perday = r * b * I_list[i - 1] * S_list[i - 1] / N E_perday_for_E = r2 * b2 * E_list[i - 1] * S_list[i - 1] / N I_perday = a * E_list[i - 1] R_perday = y * I_list[i - 1] D_perday = d * I_list[i - 1] C_list.append(C_list[i - 1] + I_perday) I_list.append(I_list[i - 1] + I_perday - R_perday - D_perday) S_list.append(S_list[i - 1] - E_perday - E_perday_for_E) R_list.append(R_list[i - 1] + R_perday) D_list.append(D_list[i - 1] + D_perday) E_list.append(E_list[i - 1] + E_perday + E_perday_for_E - I_perday)
我们建立一个 SEIR(D) 4.0 模型示例看看，设置的参数如下所示：
- 区域内的总人口
- 感染者人数
- 潜伏者人数
- 每个感染者每天接触的人数
- 是戒严后的 ，
- 每个潜伏者每天接触的人数
- 是戒严后的 ，
- 从第几天开始戒严
- 感染者将病毒传染给易感者的概率
- 潜伏者将病毒传染给易感者的概率
- 口罩防护率
- 每天感染者被治愈的人数占感染者总数的比例
- 每天感染者死亡的人数占感染者总数的比例
- 潜伏者每天转化为感染者的概率
- 模型需要模拟的总天数为
我们使用 Python 3.8 和 ReportLab 进行模拟：
上图中：The Infected 是感染者人数，Susceptibles 是易感者人数，The Recovered 是康复者人数，Death 是病死者人数，The Exposed 是潜伏者人数，Total Population 是总人口。
我们可以发现，相比于 SEIR 3.0 模型示例，在 SEIR(D) 4.0 模型示例中疫情得到了进一步的控制：
- 疾病高发期的时间向后延长了大约 20 天
- 潜伏者和感染者高峰人数下降至 2000 人以内
当然，SEIR 3.0 模型中没有病死者这一变量。
通过上述两个模型我们能够得出结论：封城（或者相对其他一些国家来说是宣布国家进入紧急状态）是行之有效的办法；如果非得出门，还请戴口罩。
至此，所有模型就介绍完了，现在我们再举几个模型示例。
模型举例
我们现在使用 SEIR 2.0 模型来看看宅在家里的重要性，设置的参数如下所示：
- 区域内的总人口
- 感染者人数
- 潜伏者人数
- 每个感染者每天接触的人数 、 和
- 每个潜伏者每天接触的人数 、 和
- 感染者将病毒传染给易感者的概率
- 潜伏者将病毒传染给易感者的概率
- 每天感染者被治愈的人数占感染者总数的比例
- 潜伏者每天转化为感染者的概率
- 模型需要模拟的总天数为
其中 和 的值会不同，分别取值 3 、10 和 50 。
通过上述三幅图我们能够明显观察到宅在家里不出门的重要性，每天大家接触的人越少，疫情就越容易得到控制。
我们再使用 SEIR 2.0 模型来看看戴口罩的重要性，设置的参数如下所示：
- 区域内的总人口
- 感染者人数
- 潜伏者人数
- 每个感染者每天接触的人数
- 每个潜伏者每天接触的人数
- 感染者将病毒传染给易感者的概率 和
- 潜伏者将病毒传染给易感者的概率 和
- 每天感染者被治愈的人数占感染者总数的比例
- 潜伏者每天转化为感染者的概率
- 模型需要模拟的总天数为
其中 和 的值会不同，分别取值 0.03 和 0.003 。
我们可以看到，当概率 和 下降到原来的十分之一时，疫情在很大程度上也能得到控制。
接下来是对新冠肺炎的简单预测。
对新冠肺炎的简单预测
Ricky 将会使用 SEIR(D) 4.0 模型对此次新冠肺炎疫情做出简单的预测。武汉市的疫情最为严峻，所以 Ricky 决定对湖北省的疫情做出预测，并且做两次预测：一次是最坏的预测，一次是为了拟合真实数据而做的预测。我们现在来推算一下 SEIR(D) 4.0 模型的那 15 个参数：
1 、区域内的总人口
根据《 2018 年湖北省国民经济和社会发展统计公报》年末全省常住人口 5917 万人，人口自然增长率为 4.54 ‰ ，假设 2019 年人口增长率为 4.5‰ ，则 2019 年 12 月湖北省总人口约为 5940 万人。
2 、感染者人数
3 、潜伏者人数
我们通过《关于新冠肺炎的一切｜回形针》的视频可以知道，第一位确诊患者于 2019 年 12 月 8 日就医：
我们知道大部分患者的潜伏期在 2 到 14 天，最长的潜伏期甚至有 24 天的。根据《钟南山领衔在顶级医学期刊发文，回应超长潜伏期 | 澎湃新闻》可知，新冠肺炎的平均潜伏期为 4 天。所以如果从第一例感染者开始计算，那么 2019 年 12 月 8 日是第一天；如果从第一例潜伏者开始计算，那么 2019 年 12 月 4 日是第一天。因为潜伏期也具有传染性，所以在建立模型时，我们从第一例潜伏者（也就是 2019 年 12 月 4 日）开始预测。
4 、每个感染者每天接触的人数
5 、 是戒严后的 ，
6 、每个潜伏者每天接触的人数
7 、 是戒严后的 ，
感染者一经确诊就会被医院隔离进行人身活动限制，因此，宜将每个感染者每天有效接触的平均人数（日接触率）参照对应于 SARS 疫情管控后的日接触率 6 。而潜伏者并没有被医院确诊，也可能并没有发病，尚可以自由活动，因此，每个潜伏者每天有效接触的平均人数（日接触率）参照对应于 SARS 疫情管控前的日接触率 12 。
而又因为当年交通便利程度相比于 2003 年 SARS 爆发时期更高，加之，本次疫情爆发正值春运时节，因此，考虑到这一流动性因素，假设日接触率会提高 10 % 。即感染者日接触率也上调为 6.6 ；潜伏者日接触率上调为 13.2 。
中国每平方公里平均人口密度为 143 人，湖北省平均人口密度为 325 人，是全国平均水平的 2.27 倍，加上湖北地区人口密度较大，自身交通便利，省会城市武汉又是一个九省通衢的中心城市，作为枢纽型节点，其患者和潜伏者的日接触率应该比全国水平和非湖北地区水平更高。因此，湖北地区的感染者日接触率也上调为全国水平的 2.3 倍，即 15.18 ；潜伏者日接触率上调为 30.36（参考自：《利用 SEIR 模型推演湖北、非湖北和全国疫情拐点 —— 2020 年突发风险系列》）。
但此次新冠肺炎采取的防疫措施（如封城）相较于 2003 年 SARS 疫情更为严厉，所以在 2020 年 1 月 23 日武汉市宣布封城后，这两个值还需要重新设置。目前很多小区都采取了封闭式管理，有的甚至采取了一户两天一次制度（即一户人家两天内只允许一个人出去一次，身边有朋友所在的小区也采取了该制度，网络上也有相关消息，如：《海口市龙华区：每 2 天允许 1 户 1 次 1 人进出 | 海南日报》）；而且很多城市都要求市民在乘坐公共交通如出租车、公交车和地铁时，必须戴口罩，拒不配合者可劝离（如：《杭州地铁公告：乘客必须戴口罩、测体温，否则劝离！| 杭州日报》）。所以我们可以将戒严后的每个潜伏者每天接触的人数 设置为 1 。
而在戒严后，感染者一旦确诊至少大部分都会住入隔离病房。医护人员除了戴 N95 口罩以外，还会佩戴护目镜和穿防护服等等（如下视频所示），所以 Ricky 决定将戒严后的每个潜伏者每天接触的人数 设置为 0.05（即每接触 20 个人中真正有效接触的人数为 1 人，那就是 ）。
8 、从第几天开始戒严
2020 年 1 月 23 日武汉市封城，而从 2019 年 12 月 4 日至 2020 年 1 月 23 日相差了 50 天。如果把 2019 年 12 月 4 日当做第 1 天的话，那么 2020 年 1月 23 日就是第 51 天（注意：不是第 50 天，是相差了 50 天）。
9 、感染者将病毒传染给易感者的概率
10 、潜伏者将病毒传染给易感者的概率
事先声明：Ricky 这边收集了由卫生应急办公室、湖北省卫生健康委员会和武汉市卫生健康委员会发布的几乎所有疫情通报，但因为这些数据经常是今天发布了到明天又会核减或核增前一天的数据，很多数据 Ricky 这边都没有做核减或核增，所以数据会存在误差，而基于这些数据做的推算那肯定也是存在一定的误差的。Ricky 这边收集到的数据如下图所示：
上图中：Confirmed 是累计确诊人数，Suspected 是现有疑似病例，Death 是累计死亡人数，Severe 是在治重症人数，Cure 是累计治愈人数。Death% 是病死率，Cure% 是治愈率。
疾病的传染概率目前还没有官方的统一数据。根据《新型冠状病毒感染的肺炎可疑暴露者和密切接触者管理方案（第二版）》对病例密切接触者的定义，利用新增确诊病例与累计追踪到密切接触者的比值作为替代。截至 2020 年 2 月 27 日，Ricky 这边推算的数据大约是 1.4750 %（可能有误）。
中国工程院院士、国家呼吸系统疾病临床医学研究中心主任、高级别专家组组长钟南山在接受新华社采访时曾表示：有些病人发展会比较慢，潜伏的带病毒者有多大的传染性，需要做一些观察及研究。对潜伏的带病毒者还是要注意。基于目前没有明显证据现实潜伏者与感染者的传染能力的差异，我们假设二者的传染能力相同。
11 、口罩防护率
从《关于新冠肺炎的一切｜回形针》的视频中我们可以了解到（如下图所示），普通医用口罩的微粒捕集效率大约是 80 % 。
12 、每天感染者被治愈的人数占感染者总数的比例
截至 2020 年 2 月 25 日，Ricky 这边推算的数据大约是 1.3211 %（可能有误）。
13 、每天感染者死亡的人数占感染者总数的比例
截至 2020 年 2 月 25 日，Ricky 这边推算的数据大约是 0.3888 %（可能有误）。
14 、潜伏者每天转化为感染者的概率
在数学意义上等同于潜伏期时长的倒数，而潜伏期我们假定为 4 天，所以 。
15 、模型需要模拟的总天数为
基于上述参数，我们推算出来的 SEIR(D) 4.0 模型如下图所示：
这个数据有点可怕，我们暂且认为这是最坏的打算。所以 Ricky 在上述 15 个参数的基础上修改了其中两个参数，让模型更贴合目前真实的数据：
- 每个感染者每天接触的人数
- 每个潜伏者每天接触的人数
这个模型的预测就非常准确了，我们将真实数据和模型预测出来的数据画在同一幅图里对比一下：
上图中，Confirmed 、Cure 和 Death 均为湖北省真实数据，其余数据为预测数据：
- The Infected ：当前感染者人数，是预测数据
- Total Infected ：累计感染者人数（对应累计确诊人数），是预测数据
- Confirmed ：累计确诊人数，是湖北省真实数据
- PCure ：累计康复者人数（累计治愈人数），是预测数据
- Cure ：累计康复者人数（累计治愈人数），是湖北省真实数据
- PDeath ：累计病死者人数（累计死亡人数），是预测数据
- Death ：累计病死者人数（累计死亡人数），是湖北省真实数据
- Current Day ：如果 2019 年 12 月 4 日是第 1 天的话，那么 2020 年 2 月 27 日就是第 86 天，目前真实数据仅采集到 2020 年 2 月 27 日
这里需要注意的是，Total Infected 是累计感染者人数，而 Confirmed 是累计确诊人数，患者经过潜伏期开始发病后（如咳嗽等）就属于感染者，而感染者去医院就诊甚至确诊都是需要一段时间的，所以 Total Infected 这条曲线会比 Confirmed 这条曲线提前上升，差了十天左右。
我们从中可以看到：
- 在第 86 天，The Infected 当前感染者人数在 35000 多一点。根据湖北省卫健委 2020 年 2 月 28 日的通报，截至 27 日 24 时仍在院治疗 32878 例，但这 32878 例是确诊人数，可能还有一些感染者是没有确诊的，如何印证这一点？我们可以看后些天湖北省的通报：截至 28 日 24 时全省新增新冠肺炎确诊病例 423 例，截至 29 日 24 时全省新增新冠肺炎确诊病例 570 例。
- PDeath 预测的累计死亡人数比 Death 真实的累计死亡人数要高一点，这个不好评价。据说中国的疑似病例还未确诊就死亡不算在累计死亡人数里，只有确诊病例死亡了才算入累计死亡人数《如果疑似病例还未确诊就死亡，这个死亡人数如何统计？| 知乎》。
- Cure 真实的累计治愈人数比 PCure 预测的累计治愈人数上升的幅度要快，也算是个好消息。
对新冠肺炎的简单预测（持续更新部分）
疫情还未完全结束，很多参数还在变动，所以关于新冠肺炎的简单预测还会不定期更新。
2020 年 03 月 03 日
目前，湖北省的治愈率已经超过 50 % ，也算是一个好消息了。但意大利、韩国和伊朗已经开始爆发，而且有从国外回国然后确诊的输入病例了。
不过我感觉韩国还好，根据 2020 年 02 月 28 日的数据来看：韩国死亡 13 例，确诊 2022 例，死亡率 0.6 % ，说明分母比较大，虽然疫情逐渐扩散，但至少大部分轻症都确诊了。而伊朗感觉有点严重，根据 2020 年 02 月 28 日的数据来看：死亡 26 例，确诊 270 例，死亡率 9.6 % ，怕是到重症才收治的？
我们再来看看为了拟合真实数据而做的预测（部分参数有改动）：
主要参考自：
发表评论?