使用J语言建模COVID-19疫情

# 用 J 语言建模 COVID-19 疫情 来源：https://datakinds.github.io/2020/03/15/modeling-the-coronavirus-outbreak-with-j *前言：我并非流行病学家，也不自称是。然而，我是一名数学专业出身的研究人员，已尽到审慎报告、确保准确无误的责任。文末附有完整代码，若您不想逐行翻看数据，可直接跳至结尾。* 如果你一直在关注新闻，想必已听说 COVID-19 大流行。鉴于美国疾病控制与预防中心（CDC）发布模型预测 62% 的美国人将受到感染（https://thehill.com/homenews/state-watch/487489-worst-case-coronavirus-models-show-massive-us-toll）¹，我开始好奇他们是如何得出这些数字的。事实上，我们究竟是如何对疾病传播得出任何结论的？ 已有大量成熟的疾病传播模型。最简单的模型是 SIR 模型（https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model）⁵，但它存在一些局限。第一个局限是，它假设一旦疾病在一个人体内走完病程，此人将永久免疫。第二个局限是，该模型假设所有感染者立即出现症状。这两条假设对 COVID-19 均不成立：该疾病的二次感染率为 14%（https://www.forbes.com/sites/brucelee/2020/03/15/can-you-get-infected-by-coronavirus-twice-how-does-covid-19-immunity-work/）⁶，且症状出现需要 10 到 14 天。 目前用于 COVID-19 大流行的模型是 SEIRS 模型（https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SEIR_model）\*,⁴,⁵。在 SEIR 模型中，每个人可能处于四种状态：“易感”、“暴露”、“传染”和“康复”（https://www.idmod.org/docs/hiv/model-seir.html）⁷。 - 一个人如果**易感**，表示他们有可能感染疾病。 - 一个人如果**暴露**，表示他们已感染或携带疾病，但**未出现**任何症状。这发生在潜伏期。 - 一个人如果**传染**，表示他们已感染疾病且正在出现症状。这时人们会去医院、需要重症监护等。 - 一个人如果**康复**，表示他们已克服疾病，并具备一定免疫力。 我们有四个可调整的变量：\\(\beta, \sigma, \gamma, \xi\\)。它们都是以通用时间单位为单位的速率。 - \\(\beta\\) 是易感者转为暴露者的速率。 - \\(\sigma\\) 是暴露者转为传染者的速率。 - \\(\gamma\\) 是传染者康复的速率。 - \\(\xi\\) 是康复者失去免疫力转为易感者的速率。 SEIRS 模型中的状态和变量。 *图片来自 https://www.idmod.org/docs/hiv/_images/SEIR-SEIRS.png，已获授权使用。版权所有 © 2019, Intellectual Ventures Management, LLC (IVM)。保留所有权利。* 这些状态和变量都有相关的微分方程。不过仍需注意几点。首先，该模型不反映因疾病导致的死亡率。因此，疾病的基础再生数（https://en.wikipedia.org/wiki/Basic_reproduction_number）⁴,⁷（即平均一人在整个传染期内会感染的人数）计算起来很简单：\\(R_0 = \frac{\beta}{\gamma}\\)。 不过，我们打算稍微复杂化一点。只看疾病数字太无聊，也不够直观！我们不分析微分方程，而是看一个实际的离散时间模拟，它使用了这四个变量。但首先，我们需要确定这些变量应该取什么值。 我们将“时间单位”设为一天。以下是关于该病毒的一些已知信息： - 我们知道 COVID-19 的 \\(R_0\\) 约为 2.28（https://www.ncbi.nlm.nih.gov/pubmed/32097725）⁸。 - 我们知道潜伏期约为 10 天左右（2-14 天，根据 CDC (https://www.cdc.gov/coronavirus/2019-ncov/symptoms-testing/symptoms.html)）。 - 我们知道一旦出现症状，症状持续约两周。⁹ - 我们知道 14%⁶ 的感染者会再次感染。 根据这些信息，我们可以计算所有需要的变量，首先从 \\(\sigma\\) 开始。每天，COVID-19 开始出现症状有一定概率 \\(\sigma\\)。每天，它*累计***未**出现症状的概率是 \\((1 - \sigma)^N\\)，其中 \\(N\\) 是天数。这是一个简单的二项分布，要使它在平均 10 天后开始出现症状，只需解 \\((1 - \sigma)^{10} = 0.5\\)。因此，\\(\sigma = 0.066967\\)。 一旦症状出现，持续约两周。我们采用相同过程：\\((1 - \gamma)^{14} = 0.5\\) 得到 \\(\gamma = 0.0483048\\)。 我们可以利用前面确定的 \\(R_0\\) 和 \\(\gamma\\) 值⁸，通过 \\(R_0 = \frac{\beta}{\gamma}\\) 计算 \\(\beta\\)。\\(2.28 = \frac{\beta}{0.0483048}\\) 得 \\(\beta = 0.110134944\\)。 剩下的就是计算 \\(\xi\\)。我不知道如何准确计算，但我们不妨假设（完全随机、合理的猜测）失去 COVID-19 免疫力需要约 30 天。\\((1 - \xi)^{30} = 0.5\\) 给出 \\(\xi = 0.02284\\)。 现在所有变量都已确定，让我们启动 J 会话，开始编写这个模拟程序！ 首先，数值定义。 `` ppl =: 10 NB. 模拟中的人数？ Beta =: 0.110134944 NB. 易感者转为暴露者的速率 Sigma =: 0.066967 NB. 暴露者转为传染者的速率 Gamma =: 0.0483048 NB. 传染者康复的速率 Xi =: 0.02284 NB. 康复者失去免疫力转为易感者的速率 R =: 2.28 NB. 一个人会感染多少人？ CFR =: 0.02 NB. 病死率 D =: 0.00155285 NB. 出现症状期间每天的死亡率 `` 注意，我们额外添加了一个参数，即病死率（https://wwwnc.cdc.gov/eid/article/26/6/20-0320_article）⁹。这是 COVID-19 病例中最终死亡的比例。CDC 建议使用 0.25% 到 3%，所以我取 2% 作为合理中间值，尤其是在美国医疗系统将不堪重负的情况下。然后，利用该值以及从症状到死亡需要 13 天的事实⁹，通过公式 \\(0.98 = (1-D)^{13}\\) 得到出现症状期间每天死亡的概率 \\(D = 0.00155285\\)。 如果我们想要模拟疾病的传播，就必须模拟社交接触。因此，让我们创建一个厄米矩阵（https://en.wikipedia.org/wiki/Hermitian_matrix）来表示人们之间的亲密度。 `` ] closeness =: 2 %~ (+ |:) ? (ppl,ppl) $ 0 0.989266 0.433774 0.587923 0.860137 0.423383 0.669393 0.296438 0.658781 0.253612 0.327771 0.433774 0.797844 0.412043 0.734478 0.693995 0.452715 0.585249 0.67874 0.591025 0.897202 0.587923 0.412043 0.916071 0.336036 0.652928 0.531987 0.541966 0.477779 0.522089 0.293315 0.860137 0.734478 0.336036 0.125376 0.438783 0.451785 0.695368 0.468632 0.508739 0.603744 0.423383 0.693995 0.652928 0.438783 0.557975 0.450049 0.722475 0.560327 0.676762 0.809449 0.669393 0.452715 0.531987 0.451785 0.450049 0.435358 0.839259 0.223557 0.835785 0.611344 0.296438 0.585249 0.541966 0.695368 0.722475 0.839259 0.380239 0.36727 0.321142 0.258312 0.658781 0.67874 0.477779 0.468632 0.560327 0.223557 0.36727 0.0571881 0.41249 0.63065 0.253612 0.591025 0.522089 0.508739 0.676762 0.835785 0.321142 0.41249 0.483952 0.657467 0.327771 0.897202 0.293315 0.603744 0.809449 0.611344 0.258312 0.63065 0.657467 0.83782 `` 这个矩阵表示第 \\(i\\) 个人和第 \\(j\\) 个人之间进行社交接触的可能性。由于自己给自己传播病毒不合理，我们去掉主对角线。 `` ] risk =: closeness - closeness * =i.ppl 0 0.433774 0.587923 0.860137 0.423383 0.669393 0.296438 0.658781 0.253612 0.327771 0.433774 0 0.412043 0.734478 0.693995 0.452715 0.585249 0.67874 0.591025 0.897202 0.587923 0.412043 0 0.336036 0.652928 0.531987 0.541966 0.477779 0.522089 0.293315 0.860137 0.734478 0.336036 0 0.438783 0.451785 0.695368 0.468632 0.508739 0.603744 0.423383 0.693995 0.652928 0.438783 0 0.450049 0.722475 0.560327 0.676762 0.809449 0.669393 0.452715 0.531987 0.451785 0.450049 0 0.839259 0.223557 0.835785 0.611344 0.296438 0.585249 0.541966 0.695368 0.722475 0.839259 0 0.36727 0.321142 0.258312 0.658781 0.67874 0.477779 0.468632 0.560327 0.223557 0.36727 0 0.41249 0.63065 0.253612 0.591025 0.522089 0.508739 0.676762 0.835785 0.321142 0.41249 0 0.657467 0.327771 0.897202 0.293315 0.603744 0.809449 0.611344 0.258312 0.63065 0.657467 0 `` 现在，让我们创建我们的人群！用 `0` 表示已死亡的人，`1` 表示易感者，`2` 表示暴露者，`3` 表示传染者，最后 `4` 表示康复者。 我们的人群可以设为 10 个易感者加一名感染者。 `` ] starting_pop =: 2 , (ppl - 1) $ 1 2 1 1 1 1 1 1 1 1 1 `` 接下来定义 `infect` 函数，给定一个人群，它使疾病在单位时间步长内推进。它将是一个一元（单参数）函数，参数是当前人群： 我们现在进入 `infect` 的定义。首先获取一个布尔向量，代表人群中谁具有传染性。 `` can_spread =. (1&< *. 4&>) y `` 将该向量应用于风险矩阵，使得每一行代表第 `j` 列中 `n` 个人与一名感染者进行社交接触的总和。 `` infectiousness =. can_spread ,./ . * risk `` 让我们看看目前的结果（为了清晰起见，添加了输出语句）。 `` infect =: 3 : 0 smoutput ] can_spread =. (1&< *. 4&>) y smoutput ] infectiousness =. can_spread ,./ . * risk ) starting_pop 2 1 1 1 1 1 1 1 1 1 risk 0 0.320427 0.628024 0.856909 0.705072 0.162919 0.684089 0.0256438 0.531406 0.88479 0.320427 0 0.653424 0.152774 0.607565 0.589595 0.232512 0.534223 0.724831 0.739889 0.628024 0.653424 0 0.485708 0.693178 0.307283 0.455467 0.381945 0.193344 0.431928 0.856909 0.152774 0.485708 0 0.748798 0.0610626 0.490139 0.803612 0.724202 0.437314 0.705072 0.607565 0.693178 0.748798 0 0.658209 0.722431 0.330002 0.355295 0.69828 0.162919 0.589595 0.307283 0.0610626 0.658209 0 0.306342 0.590838 0.779528 0.632452 0.684089 0.232512 0.455467 0.490139 0.722431 0.306342 0 0.820687 0.481649 0.315347 0.0256438 0.534223 0.381945 0.803612 0.330002 0.590838 0.820687 0 0.638455 0.561562 0.531406 0.724831 0.193344 0.724202 0.355295 0.779528 0.481649 0.638455 0 0.770651 0.88479 0.739889 0.431928 0.437314 0.69828 0.632452 0.315347 0.561562 0.770651 0 infect starting_pop 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.320427 0 0 0 0 0 0 0 0 0 0.628024 0 0 0 0 0 0 0 0 0 0.856909 0 0 0 0 0 0 0 0 0 0.705072 0 0 0 0 0 0 0 0 0 0.162919 0 0 0 0 0 0 0 0 0 0.684089 0 0 0 0 0 0 0 0 0 0.0256438 0 0 0 0 0 0 0 0 0 0.531406 0 0 0 0 0 0 0 0 0 0.88479 0 0 0 0 0 0 0 0 0 `` 我们看到 `infect starting_pop` 打印出第 1 个人具有传染性，并且他们使一些人面临很高的风险——具体来说，第 4 个人的接触等级为 `0.856909`，第 10 个人的接触等级为 `0.88479`。 在我们的模拟中，人类接触的平均水平为 `0.5`，因此我们将以此作为基线。任何高于 `0.5` 的值会增大 \\(\beta\\)，任何低于 `0.5` 的值会减小 \\(\beta\\)。这模拟了现实情况，即人们往往更容易将疾病传播给与他们关系更密切的人。 现在，让我们分别计算人群中每个人的风险。这里我们处理四种风险，每种风险可以用相应的概率向量表示。首先，将人们分为这些组： `` infect =: 3 : 0 can_spread =. (1&< *. 4&>) y susceptible =. y = 1 exposed =. y = 2 infectious =. y = 3 recovered =. y = 4 smoutput ] infectiousness =. can_spread ,./ . * risk ) `` 现在找出每个人都有风险的社交接触水平，并以某种方式对其进行限制。这有点模糊……如果一个人的所有患病朋友的亲密度评分都是 1，会发生什么？流行病学模型通常不需要在模拟人群时考虑社交亲密度，因为它们通常不模拟个体本身。因此，我将取一个人所有社交接触值的平均值，除以可能感染他们的总人数。这实际上将一个人的社交接触值转换为一个二项分布，最大值为 `1`，最小值为 `0`，众数为 `0.5`。不过，我更希望将其缩放，使其最小值约为 `-0.1`，最大值约为 `0.1`，众数为 `0`，这样社交/非社交性并非决定一个人是否生病的唯一因素。如下所示： `` 5 %~ 0.5 -~ (+/ can_spread) %~ +/ |: susceptible * infectiousness `` 函数的其余部分很简单，我们只需添加每个类别中人员的速率，然后为每个人每天计算一些随机数，如果他们不幸，就推动疾病进程。以下是完整的 `infect` 函数。 `` infect =: 3 : 0 can_spread =. (1&< *. 4&>) y susceptible =. y = 1 exposed =. y = 2 infectious =. y = 3 recovered =. y = 4 infectiousness =. can_spread ,./ . * risk susceptible_to_exposed =. (Beta * susceptible) + 5 %~ 0.5 -~ (+/ can_spread) %~ +/ |: susceptible * infectiousness exposed_to_infectious =. Sigma * exposed infectious_to_recovered =. Gamma * infectious infectious_to_dead =. D * infectious recovered_to_susceptible =. Xi * recovered 1 (I.recovered_to_susceptible>?ppl$0)} 4 (I.infectious_to_recovered>?ppl$0)} 0 (I.infectious_to_dead>?ppl$0)} 3 (I.exposed_to_infectious>?ppl$0)} 2 (I.susceptible_to_exposed>?ppl$0)} y ) `` 最后，创建一个函数来打印我们的人群。 `` infect_display =: 3 : 0 out =. infect y smoutput out out ) `` 现在终于可以模拟我们的人群了！让我们以 10 个人运行模拟 100 天。 `` starting_pop 2 2 1 1 1 1 1 1 1 1 infect_display^:100 starting_pop 2 2 1 1 2 1 2 1 1 1 2 2 1 1 2 1 3 1 1 1 2 2 1 1 2 1 3 1 1 2 2 2 2 1 2 1 3 1 1 2 3 2 2 1 2 1 3 2 1 2 3 2 2 1 2 1 3 2 2 2 3 2 2 1 2 1 3 2 2 2 3 2 2 1 2 1 3 2 2 2 4 2 2 2 2 1 3 3 2 2 4 3 2 2 2 1 3 3 2 2 4 3 2 2 2 1 3 3 2 2 4 3 2 2 2 1 4 3 2 2 4 3 2 2 2 1 4 3 2 2 4 3 2 2 2 1 4 3 2 2 4 3 2 2 2 1 4 3 2 2 4 3 2 2 2 1 4 3 2 2 4 3 2 2 2 1 4 3 2 2 4 3 2 2 2 1 4 3 2 2 4 3 2 2 2 1 4 3 2 2 4 3 2 2 2 1 4 3 2 2 4 3 3 2 2 1 4 3 2 2 4 3 3 2 2 1 4 3 2 2 4 4 3 2 2 2 4 3 2 2 4 4 3 2 2 2 4 3 2 2 4 4 3 2 2 2 4 3 2 2 4 4 3 2 3 2 4 3 2 2 4 4 3 3 3 2 4 4 2 2 4 4 3 4 3 2 4 4 2 2 4 4 1 4 3 2 4 4 3 2 4 4 1 4 3 2 4 4 3 2 4 4 2 4 3 2 4 4 3 2 4 4 2 4 3 2 4 4 3 2 4 4 2 4 3 2 4 4 3 2 4 4 2 4 3 2 4 4 3 2 4 4 2 4 0 2 4 4 3 2 4 4 2 4 0 2 4 4 3 2 4 4 2 4 0 2 4 4 3 2 4 4 2 4 0 2 4 4 4 2 4 4 2 4 0 2 4 4 4 2 4 4 2 4 0 2 4 4 4 2 4 4 2 4 0 2 4 4 4 2 4 4 2 4 0 2 4 4 4 2 4 4 3 4 0 2 4 4 4 2 4 4 3 4 0 2 4 4 4 2 4 4 3 4 0 2 4 4 4 2 4 4 4 4 0 3 4 4 4 2 4 4 4 4 0 3 4 4 4 2 4 4 4 4 0 3 4 4 4 2 4 4 4 4 0 1 4 4 4 3 4 4 4 4 0 1 4 4 4 3 4 4 4 4 0 1 4 4 4 3 4 4 4 4 0 1 4 4 4 3 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1 4 4 4 4 0 1 4 4 4 1