分类:考虑人口流动的传染病模型

背景

考虑到新冠肺炎的特殊情况——从一个城市通过人口流动传染全国并且在全国范围内实行了强制限流干预（尤其是对于重点人口流动来源城市），构建下面的粗糙指数型模型和更精确的传染病模型。先考虑一个城市的情况（其他城市当做外界），再一起考虑全国各个城市。针对这个背景，本研究主要解决：现在湖北内外的病人数量和传染病参数怎样，未来会怎样，以及，将来考虑到类似的事情发生的时候（D城市发生某种传染能力的传染病），其他各个城市应该准备什么样的资源来应对这样的风险。

模型

简单模型和再生数的计算，教科书内容

单一地区，均匀混合，不考虑治愈，考虑传播周期内不同时间点具有不同传播能力的迭代模型和微分方程模型。由于考虑了不同时间点的传播能力，

变量定义

1. [math]\displaystyle{ m\left(t\right) }[/math]：[math]\displaystyle{ t }[/math]时刻新增被感染人数
2. [math]\displaystyle{ \beta\left(t,\tau\right) }[/math]^[1]（或者定义为[math]\displaystyle{ \tilde{\beta}\left(t,\tau\right)=\beta\left(t-\tau,\tau\right) }[/math]把被感染时间单独提出来）：在[math]\displaystyle{ t-\tau }[/math]时刻被感染的病人过了[math]\displaystyle{ \tau }[/math]时间，也就是到了[math]\displaystyle{ t }[/math]时刻的时候的传染能力（传染人数）。通常进一步假设[math]\displaystyle{ t,\tau }[/math]可分离[math]\displaystyle{ \beta\left(t,\tau\right)=R\left(t\right)w\left(\tau\right) }[/math]，也就是不管在什么时间被感染，其每个时间点的相对传染能力[math]\displaystyle{ w\left(\tau\right) }[/math]一致。同时，这也意味着[math]\displaystyle{ w\left(\tau\right) }[/math]已经被归一化，也就是[math]\displaystyle{ \sum_{\tau=0}^{\infty}w\left(\tau\right)=1 }[/math]。在实际中，通常，[math]\displaystyle{ w\left(\tau\right) }[/math]在[math]\displaystyle{ \tau }[/math]很大的时候等于零，也就是传播期有限；同时，很可能[math]\displaystyle{ w\left(\tau\right) }[/math]会分成不同的阶段，例如无传染能力的时间段（无传染能力潜伏期）和有传染能力的时间段（有传染能力的潜伏期和发病期，甚至治愈之后一段时间）。假设有一个病人在生存周期内一直存活（注意本模型不考虑治愈，或者说治愈的效果已经被放到了[math]\displaystyle{ w\left(\tau\right) }[/math]里面），那么，其传染的病人数量为[math]\displaystyle{ \sum_{\tau=0}^{\infty}\beta\left(t,\tau\right)=\sum_{\tau=0}^{\infty}R\left(t\right)w\left(\tau\right)=R\left(t\right) }[/math]。

基本模型

1. [math]\displaystyle{ m\left(t\right)=\sum_{\tau=1}^{\infty}m\left(t-\tau\right)R\left(t\right)w\left(\tau\right) }[/math]

根据这个模型，我们可以得到[math]\displaystyle{ R\left(t\right)=\frac{m\left(t\right)}{\sum_{\tau=1}^{\infty}m\left(t-\tau\right)w\left(\tau\right)} }[/math]。

在这个模型下，一切疾病本身的传染能力、防控措施、治愈等因素，都被放到了[math]\displaystyle{ \beta\left(t,\tau\right)=R\left(t\right)w\left(\tau\right) }[/math]里面。在这个框架下，再生数[math]\displaystyle{ R\left(t\right) }[/math]的计算也非常简单，只需要[math]\displaystyle{ t }[/math]时刻之前的病人数量时间序列[math]\displaystyle{ m\left(t-\tau\right) }[/math]和传染期传染能力分布函数[math]\displaystyle{ w\left(\tau\right) }[/math]就可以。后者可以通过调查每一个病人在被传染的时候是传染者病程的哪一天来得到。当然，也可以简单估计出来均值和方差，然后运用典型的分布函数来近似得到这个传染能力函数。上面这个计算过程已经被编成一个计算软件包EpiEstim^[2]。

具体到本次传染病，由于湖北之内的病人数量没有准确的数据，只能通过湖北之外的数据来大概估计了。因此，在全国去除湖北的数据上做这个再生数时间序列的计算就可以。

在传染病初期，病人数量通常符合指数增长[math]\displaystyle{ m\left(t\right)=m_{0}e^{\gamma t} }[/math]。代入再生数一般公式，我们消掉了时间[math]\displaystyle{ t }[/math]，于是，在疾病初期，得到[math]\displaystyle{ R_{0}=\frac{1}{\sum_{\tau=1}^{\infty}e^{-\gamma\tau}w\left(\tau\right)} }[/math]。

如果我们假设[math]\displaystyle{ w\left(\tau\right)=\frac{1}{T} }[/math]在周期[math]\displaystyle{ T }[/math]的范围内是一个常数，于是[math]\displaystyle{ R\left(t\right)=\frac{Tm\left(t\right)}{\sum_{\tau=1}^{T}m\left(t-\tau\right)} }[/math]。

如果我们进一步结合上面两个假设，则[math]\displaystyle{ R_{0}=\frac{T}{\sum_{\tau=1}^{T}e^{-\gamma\tau}}=\frac{T\left(1-e^{-\gamma}\right)}{1-e^{-\gamma\left(T+1\right)}} }[/math]。在[math]\displaystyle{ \gamma\ll 1 }[/math]的时候，得到近似解[math]\displaystyle{ R_{0}=T\gamma }[/math]。和通常的最简单的估计公式一致。

考虑输入性病人的模型

1. 特定城市输入病人（称为零代病人）的一代感染者数量[math]\displaystyle{ m_{D}\left(t\right)=\sum_{\tau=1}^{\infty}D\left(t-\tau\right)R_{D}\left(t\right)\sum_{\xi=0}^{\infty}\rho_{D}\left(\xi\right)w\left(\tau+\xi\right) }[/math]
2. 其他城市输入病人（称为零代病人）的一代感染者数量[math]\displaystyle{ m_{O}\left(t\right)=\sum_{\tau=1}^{\infty}O\left(t-\tau\right)R_{O}\left(t\right)\sum_{\xi=0}^{\infty}\rho_{O}\left(\xi\right)w\left(\tau+\xi\right) }[/math]
3. 本地城市二代以及二代以上以上感染者数量[math]\displaystyle{ m_{C}\left(t\right)=\sum_{\tau=1}^{\infty}m\left(t-\tau\right)R_{C}\left(t\right)w\left(\tau\right) }[/math]
4. [math]\displaystyle{ m\left(t\right)=m_{D}\left(t\right)+m_{O}\left(t\right)+m_{C}\left(t\right) }[/math]

其中，特定城市的输入性病人和其他城市的输入性病人可能很难得到其处于传染周期哪一天的数据。因此，加上一个随机数。这个随机数的简单处理方法就是用例如传染周期的平均值，或者说有行动能力（例如未发病）的传染期的平均值[math]\displaystyle{ T_{c}=\sum_{\tau}\tau w\left(\tau\right) }[/math]。

为了简化记号，我们定义[math]\displaystyle{ \omega_{D}\left(\tau\right)=\sum_{\xi=0}^{\infty}\rho_{D}\left(\xi\right)w\left(\tau+\xi\right) }[/math]，[math]\displaystyle{ \omega_{O}\left(\tau\right)=\sum_{\xi=0}^{\infty}\rho_{O}\left(\xi\right)w\left(\tau+\xi\right) }[/math]。

如果有了[math]\displaystyle{ R_{D}\left(t\right),R_{O}\left(t\right),R_{C}\left(t\right) }[/math]以及初始条件（任意长度的[math]\displaystyle{ m_{D}\left(t\right), m_{O}\left(t\right), m_{C}\left(t\right) }[/math]的历史数据）以及外生变量[math]\displaystyle{ D\left(t\right), O\left(t\right) }[/math]。我们就可以拟合出来历史上的和外来的[math]\displaystyle{ m_{D}\left(t\right), m_{O}\left(t\right), m_{C}\left(t\right) }[/math]。对于未来，可以把有效再生数的已知部分的末尾的平均值当做未来的有效再生数。当然，假设[math]\displaystyle{ w\left(\tau\right) }[/math]完全已知。

于是，相应的再生数为，

1. [math]\displaystyle{ R_{D}\left(t\right)=\frac{m_{D}\left(t\right)}{\sum_{\tau=1}^{\infty}D\left(t-\tau\right)w\left(\tau+T_{c}\right)} }[/math]
2. [math]\displaystyle{ R_{O}\left(t\right)=\frac{m_{O}\left(t\right)}{\sum_{\tau=1}^{\infty}O\left(t-\tau\right)w\left(\tau+T_{c}\right)} }[/math]
3. [math]\displaystyle{ R_{C}\left(t\right)=\frac{m_{C}\left(t\right)}{\sum_{\tau=1}^{\infty}\left[m_{D}\left(t-\tau\right)+m_{O}\left(t-\tau\right)+m_{C}\left(t-\tau\right)\right]w\left(\tau\right)} }[/math]

如果我们需要从实际数据中获得[math]\displaystyle{ R_{D}\left(t\right), R_{O}\left(t\right), R_{C}\left(t\right) }[/math]的话，我们需要一段时间的[math]\displaystyle{ m_{D}\left(t\right), m_{O}\left(t\right), m_{C}\left(t\right) }[/math]实际数据，以及外生变量[math]\displaystyle{ D\left(t\right), O\left(t\right) }[/math]。当然，假设[math]\displaystyle{ w\left(\tau\right) }[/math]完全已知。

如果我们只关心[math]\displaystyle{ R_{C}\left(t\right) }[/math]的话，我们只要从所有的确诊病人中剔除所有输入性病人，也就是得到[math]\displaystyle{ m_{D}\left(t\right)+m_{O}\left(t\right)+m_{C}\left(t\right) }[/math]，然后，有[math]\displaystyle{ m_{C}\left(t\right) }[/math]的数据（本地二代以及二代以上以上感染者数量）就可以计算了。

输入性病人模型的简化情况

插入一个以上方程的简化情况，假设所有的有效再生数都一样[math]\displaystyle{ R_{D}\left(t\right)=R_{O}\left(t\right)=R_{C}\left(t\right)=R\left(t\right) }[/math]，并且忽略[math]\displaystyle{ \xi }[/math]则，把三个方程加起来，可以得到，[math]\displaystyle{ m\left(t\right)=\sum_{\tau=1}^{\infty}\left(D\left(t-\tau\right)+O\left(t-\tau\right)+m\left(t\right)\right)R\left(t\right)w\left(\tau\right) }[/math]，反过来有效再生数的计算为[math]\displaystyle{ R\left(t\right)=\frac{m\left(t\right)}{\sum_{\tau=1}^{\infty}\left(D\left(t-\tau\right)+O\left(t-\tau\right)+m\left(t\right)\right)w\left(\tau\right)} }[/math]。这实际上是EpiEstim2软件^[3]包完成的对EpiEstim的拓展。

其中[math]\displaystyle{ \left(D\left(t\right)+O\left(t\right) }[/math]被称为输入病例，[math]\displaystyle{ m\left(t\right) }[/math]被称为本地病例。按照我们的语言，前者是零代病例，后者为一代以及一代以上病例。只有当前面的假设成立的时候，“输入和本地的有效再生数一样”、“输入病人都在病程零点”，这个公式才成立。更一般的情况（由于防控措施导致几种有效再生数不一样、需要考虑输入病人的实际病程），需要区分零代、一代、二代和二代以上。

插入完毕。

考虑输入性病人和被收治概率的模型

由于实际上的病人数量是不客观测量，可观测量往往是被收治的病人数量，我们引入下面的被收治概率的概念。

1. 特定城市输入病人的传染、收治
   1. 在[math]\displaystyle{ t }[/math]时刻来到所关心的城市的来自于特定城市的输入病人（称为零代病人）的数量，[math]\displaystyle{ D\left(t,0\right) }[/math]。假设当天不会被收治。如果当天被收治则考虑直接减少[math]\displaystyle{ D\left(t,0\right) }[/math]，把其中的一部分直接算作收治的输入病人里面
   2. 在[math]\displaystyle{ t }[/math]时刻来到所关心的城市的来自于特定城市的输入病人（称为零代病人）过了[math]\displaystyle{ \tau }[/math]时间以后剩下的数量，[math]\displaystyle{ D\left(t,\tau\right)=D\left(t,0\right)\left(1-\sum_{\hat{\tau}=1}^{\tau}O_{D}\left(\hat{\tau}\right)\right) }[/math]，其中[math]\displaystyle{ O_{D}\left(\tau\right) }[/math]是一个输入病人在输入之后的第[math]\displaystyle{ \tau }[/math]天被收治的概率。注意，在这里其实引入了一个近似，[math]\displaystyle{ O_{D}\left(\tau\right) }[/math]本来应该是[math]\displaystyle{ O_{D}\left(t,\tau\right) }[/math]同时依赖于输入病人进入所关心的城市的时间[math]\displaystyle{ t }[/math]和进入以后的时间[math]\displaystyle{ \tau }[/math]
   3. 在[math]\displaystyle{ t }[/math]时刻被收治的特定城市输入病人（称为零代病人）数量，[math]\displaystyle{ X_{D}\left(t\right)=\sum_{\tau=1}^{\infty} X_{D}\left(t, \tau\right)=\sum_{\tau=1}^{\infty}D\left(t-\tau,0\right)O_{D}\left(\tau\right) }[/math]
   4. 在[math]\displaystyle{ t }[/math]时刻被特定城市输入病人（称为零代病人）感染的一代感染者数量，[math]\displaystyle{ m_{D}\left(t,0\right)=\sum_{\tau}D\left(t-\tau, \tau\right)R_{D}\left(t\right)\omega_{D}\left(\tau\right) }[/math]
   5. 在[math]\displaystyle{ t }[/math]时刻被特定城市输入病人感染的一代病人过了[math]\displaystyle{ \tau }[/math]时间以后剩下的数量，[math]\displaystyle{ m_{D}\left(t,\tau\right)=m_{D}\left(t,0\right)\left(1-\sum_{\hat{\tau}=1}^{\tau}O_{m_D}\left(\hat{\tau}\right)\right) }[/math]，其中[math]\displaystyle{ O_{m_D}\left(\tau\right) }[/math]是一个一代病人在被感染之后的第[math]\displaystyle{ \tau }[/math]天被收治的概率
   6. 在[math]\displaystyle{ t }[/math]时刻被收治的一代感染者数量[math]\displaystyle{ X_{m_D}\left(t\right)=\sum_{\tau}X_{m_D}\left(t,\tau\right) = \sum_{\tau}m_{D}\left(t-\tau,\tau\right)O_{m_D}\left(\tau\right) }[/math]
2. 其他城市输入病人的传染、收治
   1. 在[math]\displaystyle{ t }[/math]时刻来到所关心的城市的来自于其他城市的输入病人（称为零代病人）的数量，[math]\displaystyle{ O\left(t,0\right) }[/math]。假设当天不会被收治。
   2. 在[math]\displaystyle{ t }[/math]时刻来到所关心的城市的来自于其他城市的输入病人（称为零代病人）过了[math]\displaystyle{ \tau }[/math]时间以后剩下的数量，[math]\displaystyle{ O\left(t,\tau\right)=O\left(t,0\right)\left(1-\sum_{\hat{\tau}=1}^{\tau}O_{O}\left(\hat{\tau}\right)\right) }[/math]，其中[math]\displaystyle{ O_{O}\left(\tau\right) }[/math]是一个其他城市的输入病人在输入之后的第[math]\displaystyle{ \tau }[/math]天被收治的概率
   3. 在[math]\displaystyle{ t }[/math]时刻被收治的其他城市输入病人（称为零代病人）数量，[math]\displaystyle{ X_{O}\left(t\right)=\sum_{\tau=1}^{\infty} X_{O}\left(t, \tau\right)=\sum_{\tau=1}^{\infty}O\left(t-\tau,0\right)O_{O}\left(\tau\right) }[/math]
   4. 在[math]\displaystyle{ t }[/math]时刻被其他城市输入病人（称为零代病人）感染的一代感染者数量，[math]\displaystyle{ m_{O}\left(t,0\right)=\sum_{\tau}O\left(t-\tau, \tau\right)R_{O}\left(t\right)\omega_{O}\left(\tau\right) }[/math]
   5. 在[math]\displaystyle{ t }[/math]时刻被其他城市输入病人感染的一代病人过了[math]\displaystyle{ \tau }[/math]时间以后剩下的数量，[math]\displaystyle{ m_{O}\left(t,\tau\right)=m_{O}\left(t,0\right)\left(1-\sum_{\hat{\tau}=1}^{\tau}O_{m_O}\left(\hat{\tau}\right)\right) }[/math]，其中[math]\displaystyle{ O_{m_O}\left(\tau\right) }[/math]是一个一代病人在被感染之后的第[math]\displaystyle{ \tau }[/math]天被收治的概率
   6. 在[math]\displaystyle{ t }[/math]时刻被收治的一代感染者数量[math]\displaystyle{ X_{m_O}\left(t\right)=\sum_{\tau}X_{m_O}\left(t,\tau\right) = \sum_{\tau}m_{O}\left(t-\tau,\tau\right)O_{m_O}\left(\tau\right) }[/math]
3. 本地城市二代以及二代以上病人的传染、收治
   1. 本地城市二代以及二代以上以上感染者数量[math]\displaystyle{ m_{C}\left(t\right)=\sum_{\tau=1}^{\infty}\left[m_{D}\left(t-\tau,\tau\right)+m_{O}\left(t-\tau,\tau\right)+m_{C}\left(t-\tau,\tau\right)\right]R_{C}\left(t\right)w\left(\tau\right) }[/math]
   2. 在[math]\displaystyle{ t }[/math]时刻被本地感染的病人过了[math]\displaystyle{ \tau }[/math]时间以后剩下的数量，[math]\displaystyle{ m_{c}\left(t,\tau\right)=m_{c}\left(t,0\right)\left(1-\sum_{\hat{\tau}=1}^{\tau}O_{m_{C}}\left(\hat{\tau}\right)\right) }[/math]，其中[math]\displaystyle{ O_{m_{C}}\left(\tau\right) }[/math]是一个本地二代以及二代以上病人在输入之后的第[math]\displaystyle{ \tau }[/math]天被收治的概率
   3. 在[math]\displaystyle{ t }[/math]时刻被收治的本地二代以及二代以上病人数量，[math]\displaystyle{ X_{m_{C}}\left(t\right)=\sum_{\tau=1}^{\infty} X_{m_{C}}\left(t, \tau\right)=\sum_{\tau=1}^{\infty}m_{c}\left(t-\tau,0\right)O_{m_{C}}\left(\tau\right) }[/math]

如何从实际数据得到收治概率？

如果[math]\displaystyle{ D\left(t,\tau\right) }[/math]是可观测的，则计算[math]\displaystyle{ 1-\frac{D\left(t,\tau\right)}{D\left(t,0\right)}=\sum_{\hat{\tau}=1}^{\tau}O_{D}\left(\hat{\tau}\right) }[/math]，于是，[math]\displaystyle{ O_{D}\left(\hat{\tau}\right) = \frac{D\left(t,\tau-1\right)-D\left(t,\tau\right)}{D\left(t,0\right)} }[/math]。当然，这个算出来的收治概率含有输入病人进入的时间[math]\displaystyle{ t }[/math]，如果进一步想得到一个不依赖于时间[math]\displaystyle{ t }[/math]，可以对所有的时间[math]\displaystyle{ t }[/math]做一个平均，也就是[math]\displaystyle{ O_{D}\left(\hat{\tau}\right) = \left\langle \frac{D\left(t,\tau-1\right)-D\left(t,\tau\right)}{D\left(t,0\right)} \right\rangle_{t} }[/math]。但是，[math]\displaystyle{ D\left(t,\tau\right) }[/math]可能不可观测。

我们假设[math]\displaystyle{ X_{D}\left(t, \tau\right) }[/math]（[math]\displaystyle{ =D\left(t-\tau,0\right)O_{D}\left(\tau\right) }[/math]，其含义是每天被收治的输入性病人中，多少是已经来到这个城市[math]\displaystyle{ \tau }[/math]天的。这个数据原则上是可以获得的）是可观测的。于是，[math]\displaystyle{ \frac{X_{D}\left(t, \tau\right)}{X_{D}\left(t-1, \tau-1\right)}=\frac{O_{D}\left(\tau\right)}{O_{D}\left(\tau-1\right)} }[/math]。得到这个比例之后，运用归一化条件[math]\displaystyle{ \sum_{\tau}O_{D}\left(\tau\right)=1 }[/math]自然可以得到[math]\displaystyle{ O_{D}\left(\tau\right) }[/math]。同样地，这个算出来的结果实际上依赖于[math]\displaystyle{ t }[/math]，可以再做一次平均来去掉这个时间[math]\displaystyle{ t }[/math]。注意平均的时候，存在多种方式，例如，[math]\displaystyle{ \left\langle\frac{X_{D}\left(t, \tau\right)}{X_{D}\left(t-1, \tau-1\right)}\right\rangle_{t} }[/math]和[math]\displaystyle{ \frac{\left\langle X_{D}\left(t, \tau\right)\right\rangle_{t}}{\left\langle X_{D}\left(t-1, \tau-1\right)\right\rangle_{t}} }[/math]，先求比例在平均和先加总（平均）再求比例。

其他的收治概率函数也可以从可观测数据也就是收治病人数据中得到。

如何从实际数据得到有效再生数？

得到[math]\displaystyle{ O_{D}\left(\tau\right) }[/math]之后，利用和标准模型一样的计算过程，就可以得到[math]\displaystyle{ R_{D}\left(t\right) }[/math]。

新冠肺炎相关计算

研究目标：大概估计一下现在的病人数量，尤其是湖北内部的病人数量。推测一下到什么时候大概疫情会结束或者说大概什么时候可以采取比现在经济代价更小的防控措施。

研究思路：考虑到湖北的历史病人数据由于检测速度和收治容量的限制等原因很难做到统计全面，我们，

1. 采用湖北之外的地区在封城之前从武汉输入的病人数量和总的武汉输入人口的比值来估计湖北内外的封城之前的病人数量；
2. 采用湖北之外的地区的本地病人有效再生数时间序列来当做湖北之内的有效再生数时间序列，然后结合第一部分的结果来推断湖北之内的病人数量的历史和未来（推测未来的时候，我们假设有效再生数已经稳定不再发生变化）；
3. 对于湖北之外的城市，我们区分输入病人有效再生数和本地病人有效再生数，来推断病人数量的历史和未来（推测未来的时候，我们假设有效再生数已经稳定不再发生变化），同时，我们讨论了对一定时间开放春运返程以及不同的返程人员防控方式下的未来趋势；

需要做的计算：

1. 全国（湖北之外）的以及个别城市的从输入性和非输入性的时间序列，估计出来的武汉封城当时总病人数量，以及按照交通流量算出来的封城当时湖北内外其他城市的病人数量
2. 按照全国湖北之外非输入性病人的时间序列（也就是分开零代、一代以及一代以上两类，然后用第二类病人数量时间序列）和[math]\displaystyle{ w(\tau) }[/math]算出来的本地病人有效再生数时间序列[math]\displaystyle{ R\left(t\right) }[/math]
   1. 如果能够分开零代、一代、二代以及二代以上，按照全国湖北之外的各类型病人数量算出来输入病人和本地病人的有效再生数时间序列[math]\displaystyle{ R_{D}\left(t\right),R_{C}\left(t\right),R_{O}\left(t\right) }[/math]。最后那个其他城市输入的有效再生数可以暂时不考虑
   2. 把这几个时间序列的趋势和来自于个案统计的平均每一个病人感染的人数的时间序列做一个对比，具体数字可以对不上，但是趋势图应该差不多一致
3. 按照全国湖北之外的本地有效再生数时间序列估计湖北病人数量的历史和未来。在这里湖北的有效再生数可以采用某种分段函数来假设。例如，从封城之后到某个时间点采用无干预再生数[math]\displaystyle{ R_{0} }[/math]，之后采用湖北之外的有效再生数。对于未来，采用有效再生数时间序列有数据的部分的最新的值（或者某个窗口的平均）。
4. 按照本地有效再生数和输入有效再生数估计其他城市的病人数量的历史和未来（在这里，可以通过对每一个城市单独重复上面的“湖北之外的全国“的计算）
5. 考虑不同城市之间的流动——本身通过交通流量来估计，来计算各个城市的未来（在这里，其他城市输入的人口可以用不同的防控方式，例如假设他们的再生数来自于湖北的一样，或者假设和本地一样，或者介于两者之间）
6. 顺便，各个城市的有效再生数的对比是否要展示？
7. 来自个案数据的传染病参数统计量的时间序列：收治之前的时间、潜伏期时间、平均传染人数等

封城带来的输入病人的特殊处理

首先，我们根据被收治的情况，得到所有的输入病人数量——也就是把每天发现输入病人加起来，记做[math]\displaystyle{ D_{0} }[/math]。于是，我们有第零天（[math]\displaystyle{ t=0 }[/math]，封城之前的一天，我们把封城当做第一天[math]\displaystyle{ t=1 }[/math]。如果需要平移时间零点，下面的公式都要做相应的变换）的输入病人数量为[math]\displaystyle{ D\left(t=0\right)=D_{0} }[/math]。在这里，我们假设所有的输入病人都是封城之前那一天进来的，并且收治这些输入病人只发生在封城当天和之后。同时，我们假设这些输入进来的病人都处于传染期的第零天（[math]\displaystyle{ \tau=0 }[/math]，如果需要可以改成一个[math]\displaystyle{ \tau }[/math]的[math]\displaystyle{ \left[0,\tau_{Max}\right] }[/math]之间的均匀分布）。然后，根据每一天被收治的输入病人的时间序列[math]\displaystyle{ \tilde{D}_{\tau} }[/math]得到输入病人数量[math]\displaystyle{ D\left(t,\hat{t}\right)=\delta_{t,0}\left(D_{0}-\sum_{1\leq \tau \leq \hat{t}}\tilde{D}_{\tau}\right) }[/math]。其中[math]\displaystyle{ \delta_{t,0} }[/math]表示输入这件事只发生在[math]\displaystyle{ t=0 }[/math]，[math]\displaystyle{ \left(D_{0}-\sum_{1\leq \tau \leq \hat{t}}\tilde{D}_{\tau}\right) }[/math]表示这些输入病人随着[math]\displaystyle{ \hat{t} }[/math]在减少（由于不断地被收治）。

注意，[math]\displaystyle{ D_{0}=\sum_{1\leq \tau \lt \infty}\tilde{D}_{\tau} }[/math]，因此，[math]\displaystyle{ D\left(t,\hat{t}\right)=\delta_{t,0}\sum_{\hat{t}\leq \tau \lt \infty} \tilde{D}_{\tau} }[/math]。

其次，有了这个输入病人的时间序列，我们代入到上面的一般公式，[math]\displaystyle{ R_{D}\left(t\right)=\frac{m_{D}\left(t\right)}{\sum_{\tau}D\left(t-\tau\right)w\left(\tau\right)}=\frac{m_{D}\left(t\right)}{w\left(t\right)\sum_{t\leq \tau \lt \infty}\tilde{D}_{\tau}} }[/math]。注意，如果时间起点换了的话，这个公式的分母部分需要重新调整。其余[math]\displaystyle{ R_{C}\left(t\right)=\frac{m_{C}\left(t\right)}{\sum_{\tau}\left(m_{D}\left(t-\tau\right)+m_{C}\left(t-\tau\right)\right)w\left(\tau\right)} }[/math]等不牵涉到输入病人的公式就不再重复了。

再次，假设我们的数据不支持我们分开[math]\displaystyle{ \left\{D\left(t\right),m_{D}\left(t\right),m_{C}\left(t\right)\right\} }[/math]（零代、一代、二代以及以上），仅仅支持我们分开[math]\displaystyle{ \left\{D\left(t\right),m\left(t\right)=m_{D}\left(t\right)+m_{C}\left(t\right)\right\} }[/math]（零代、一代以及以上），则

1. 只有在[math]\displaystyle{ m_{D}\gg m_{C} }[/math]的情况下，才能用来求[math]\displaystyle{ R_{D} }[/math]，[math]\displaystyle{ R_{D}\left(t\right)\approx \frac{m\left(t\right)}{w\left(t\right)\sum_{t\leq \tau \lt \infty}\tilde{D}_{\tau}} }[/math]；
2. 只有在[math]\displaystyle{ m_{D}\ll m_{C} }[/math]的情况下，才能用来求[math]\displaystyle{ R_{C} }[/math]，[math]\displaystyle{ R_{C}\left(t\right) \approx \frac{m\left(t\right)}{\sum_{\tau}m\left(t-\tau\right)w\left(\tau\right)} }[/math]。

当然，可以不管到底什么情况，先算出来，然后用一小部分能够拿到的数据来估计一下这个比例[math]\displaystyle{ \frac{m_{D}}{m_{c}} }[/math]。

这样算出来的[math]\displaystyle{ R_{D}\left(t\right), R_{C}\left(t\right) }[/math]，用于预测未来和拟合历史，当外生变量[math]\displaystyle{ \tilde{D}_{\tau} }[/math]完全已知的时候是比较准确的。

最后，还有一种处理方式，把输入病人数量完全当做常数[math]\displaystyle{ D\left(t\right)=\delta_{t,0}D_{0}=\delta_{t,0}\sum_{1\leq \tau\lt \infty} \tilde{D}_{\tau} }[/math]，然后不考虑每天实际剩余的输入病人数量的时间序列，反而把这个时间序列的效果当做[math]\displaystyle{ R_{D}\left(t\right) }[/math]的一部分。这样做的好处是，我们不再需要外生变量[math]\displaystyle{ \tilde{D}_{\tau} }[/math]，而仅仅需要[math]\displaystyle{ D_{0} }[/math]，而且这个[math]\displaystyle{ D_{0} }[/math]可能可以通过其他方式，例如考虑交通流OD计算出来。

到底采用哪一种处理方式，这要看后期的拟合预测任务提供的已知条件，以及对比拟合的结果。

微分方程模型（反应扩散方程）

单个城市的反应方程，考虑武汉输入和其他城市输入

1. [math]\displaystyle{ S\left(t\right) }[/math]：一个城市中的易感者人数。可以分为[math]\displaystyle{ S_{C}\left(t\right) }[/math]（本城市中的易感人数），[math]\displaystyle{ S_{D}\left(t\right) }[/math]（来自于特定城市的在本市中的易感人数），[math]\displaystyle{ S_{O}\left(t\right) }[/math]（来自于其他本城市中的在本市的易感人数）。[math]\displaystyle{ S\left(t\right)=S_{C}\left(t\right)+S_{D}\left(t\right)+S_{O}\left(t\right) }[/math]。
2. [math]\displaystyle{ L\left(t\right) }[/math]：一个城市中的潜伏者人数。可以分为[math]\displaystyle{ L_{C}\left(t\right) }[/math]（本城市中的潜伏者人数），[math]\displaystyle{ L_{D}\left(t\right) }[/math]（来自于特定城市的在本市中的潜伏者人数），[math]\displaystyle{ L_{O}\left(t\right) }[/math]（来自于其他本城市中的在本市的潜伏者人数）。[math]\displaystyle{ L\left(t\right)=L_{C}\left(t\right)+L_{D}\left(t\right)+L_{O}\left(t\right) }[/math]。
3. [math]\displaystyle{ I\left(t\right) }[/math]：一个城市中的感染者人数。可以分为[math]\displaystyle{ I_{C}\left(t\right) }[/math]（本城市中的感染者人数），[math]\displaystyle{ I_{D}\left(t\right) }[/math]（来自于特定城市的在本市中的感染者人数），[math]\displaystyle{ I_{O}\left(t\right) }[/math]（来自于其他本城市中的在本市的感染者人数）。[math]\displaystyle{ I\left(t\right)=I_{C}\left(t\right)+I_{D}\left(t\right)+I_{O}\left(t\right) }[/math]。
4. [math]\displaystyle{ X\left(t\right) }[/math]：一个城市中的被收治者人数。假设收治之后不管康复还是死亡，都不再能得病和传染。如果情况变了，将来还要在分开收治之后的康复和死亡。[math]\displaystyle{ X\left(t\right)=X_{C}\left(t\right)+X_{D}\left(t\right)+X_{O}\left(t\right) }[/math]。如果这些被收治者中如果可以区分[math]\displaystyle{ X_{CC}\left(t\right) }[/math]（本市通过接触本市人得病）、[math]\displaystyle{ X_{CD}\left(t\right) }[/math]（本市通过接触特定输入城市的输入病人得病）、[math]\displaystyle{ X_{CO}\left(t\right) }[/math]（本市通过接触其他城市的输入病人得病）、[math]\displaystyle{ X_{DC}\left(t\right) }[/math]（特定输入城市的人通过接触本市的病人得病）、[math]\displaystyle{ X_{OC}\left(t\right) }[/math]（其他城市的人通过接触本市人得病）、[math]\displaystyle{ X_{OD}\left(t\right) }[/math]（其他城市的人通过接触特定城市的病人得病）、[math]\displaystyle{ X_{DO}\left(t\right) }[/math]（特定城市的人通过接触其他城市的病人得病）、[math]\displaystyle{ X_{DC}\left(t\right) }[/math]（特定城市的人通过接触本市的病人得病）、[math]\displaystyle{ X_{DD}\left(t\right) }[/math]（特定城市的人通过特定城市的病人得病）、[math]\displaystyle{ X_{OO}\left(t\right) }[/math]（其他城市的人通过其他城市的病人得病），就更好了。注意，仅仅这里的[math]\displaystyle{ X\left(t\right) }[/math]是可以观测的。其他量目前的传染病控制方式下不能直接观测。
5. [math]\displaystyle{ Y\left(t\right) }[/math]：一个城市中的未被收治者中的恢复者和死亡人数。假设恢复和死亡之后都不再得病和传染。如果死亡当时由于尸体处理等问题会传染更多的人，则将来还要分开未收治的康复和死亡。可以分为[math]\displaystyle{ Y_{C}\left(t\right) }[/math]（本城市中的未被收治者中的恢复者和死亡人数），[math]\displaystyle{ Y_{D}\left(t\right) }[/math]（来自于特定城市的在本市中的未被收治者中的恢复者和死亡人数），[math]\displaystyle{ Y_{O}\left(t\right) }[/math]（来自于其他本城市中的在本市的未被收治者中的恢复者和死亡人数）。[math]\displaystyle{ Y\left(t\right)=Y_{C}\left(t\right)+Y_{D}\left(t\right)+Y_{O}\left(t\right) }[/math]。
6. [math]\displaystyle{ \xi_{\mu}\left(t\right) }[/math]表示从特定城市之外的其他城市输入到当地的总人口。在一定条件下，可以先不考虑这个其他城市的输入人口。

缺陷：没有直接考虑带有目标性的隔离政策（而是通过设定另一个感染率来粗略描述），例如凡是能够查到的接触过确诊病人的都隔离等措施。模型可以以后再改进。不考虑迁移到其他城市。

微分方程：

定义简便记号：

1. [math]\displaystyle{ N\left(t\right)=S\left(t\right)+L\left(t\right)+I\left(t\right) }[/math]
2. [math]\displaystyle{ S\left(t\right)=S_{C}\left(t\right)+S_{D}\left(t\right)+S_{O}\left(t\right) }[/math]
3. [math]\displaystyle{ L\left(t\right)=L_{C}\left(t\right)+L_{D}\left(t\right)+L_{O}\left(t\right) }[/math]，[math]\displaystyle{ I\left(t\right)=I_{C}\left(t\right)+I_{D}\left(t\right)+I_{O}\left(t\right) }[/math]
4. [math]\displaystyle{ X\left(t\right)=X_{C}\left(t\right)+X_{D}\left(t\right)+X_{O}\left(t\right) }[/math]，[math]\displaystyle{ I\left(t\right)=Y_{C}\left(t\right)+Y_{D}\left(t\right)+Y_{O}\left(t\right) }[/math]。后面的几项可以用矢量记号表示，例如[math]\displaystyle{ S\left(t\right)=\sum_{\mu=C,D,O}S_{\mu}\left(t\right) }[/math]。

演化方程：

1. [math]\displaystyle{ \frac{d}{dt}S_{C}\left(t\right) = -\beta^{CC}\frac{S_{C}\left(t\right)I_{C}\left(t\right)}{N\left(t\right)} -\beta^{CD}\frac{S_{C}\left(t\right)I_{D}\left(t\right)}{N\left(t\right)} -\beta^{CO}\frac{S_{C}\left(t\right)I_{O}\left(t\right)}{N\left(t\right)}-\alpha^{CC}\frac{S_{C}\left(t\right)L_{C}\left(t\right)}{N\left(t\right)}-\alpha^{CD}\frac{S_{C}\left(t\right)L_{D}\left(t\right)}{N\left(t\right)}-\alpha^{CO}\frac{S_{C}\left(t\right)L_{O}\left(t\right)}{N\left(t\right)} }[/math]，用矢量几号表示为，[math]\displaystyle{ \frac{d}{dt}S_{\mu}\left(t\right) = -\sum_{\nu}\beta^{\mu\nu}\frac{S_{\mu}\left(t\right)I_{\nu}\left(t\right)}{N\left(t\right)}-\sum_{\nu} \alpha^{\mu\nu}\frac{S_{\mu}\left(t\right)L_{\nu}\left(t\right)}{N\left(t\right)} + s_{\mu}\xi_{\mu}\left(t\right)\delta{\mu O} }[/math]。下面的方程就用矢量形式写，而不是每一项都写出来了。其中多出来的[math]\displaystyle{ s_{\mu}\xi_{\mu}\left(t\right)\delta{\mu O} }[/math]表示其他城市输入到当地的易感人群数量：[math]\displaystyle{ s_{\mu} }[/math]输入源中易感人群比例，[math]\displaystyle{ \xi_{\mu}\left(t\right) }[/math]输入总人速率。当[math]\displaystyle{ \mu=C }[/math]这一项为零。当[math]\displaystyle{ \mu=D }[/math]这一项也为零（特定城市已经封城，不再直接输出到当地）。当[math]\displaystyle{ \mu=O }[/math]这一项是一个已知时间序列（外生变量）。
2. [math]\displaystyle{ \frac{d}{dt}L_{\mu}\left(t\right) = \sum_{\nu}\beta^{\mu\nu}\frac{S_{\mu}\left(t\right)I_{\nu}\left(t\right)}{N\left(t\right)} +\sum_{\nu} \alpha^{\mu\nu}\frac{S_{\mu}\left(t\right)L_{\nu}\left(t\right)}{N\left(t\right)} -\gamma^{\mu}L_{\mu}\left(t\right)+l_{\mu}\xi_{\mu}\left(t\right)\delta{\mu O} }[/math]。其中多出来的[math]\displaystyle{ l_{\mu}\xi_{\mu}\left(t\right)\delta{\mu O} }[/math]表示其他城市输入到当地的潜伏者数量：[math]\displaystyle{ l_{\mu} }[/math]输入源中潜伏者比例。
3. [math]\displaystyle{ \frac{d}{dt}I_{\mu}\left(t\right) = \gamma^{\mu}L_{\mu}\left(t\right)-\eta^{\mu}_{X}I_{\mu}\left(t\right)-\eta^{\mu}_{Y}I_{\mu}\left(t\right)+i_{\mu}\xi_{\mu}\left(t\right)\delta{\mu O} }[/math]。其中多出来的[math]\displaystyle{ i_{\mu}\xi_{\mu}\left(t\right)\delta{\mu O} }[/math]表示其他城市输入到当地的发病人数量：[math]\displaystyle{ s_{\mu} }[/math]输入源中发病人比例。
4. [math]\displaystyle{ \frac{d}{dt}X_{\mu}\left(t\right) = \eta^{\mu}_{X}I_{\mu}\left(t\right) }[/math]，可观测量，本地病人数量 [math]\displaystyle{ X_{C}\left(t\right) }[/math]，特定城市输入病人数量 [math]\displaystyle{ X_{D}\left(t\right) }[/math]，其他城市输入病人数量 [math]\displaystyle{ X_{O}\left(t\right) }[/math]。
5. [math]\displaystyle{ \frac{d}{dt}Y_{\mu}\left(t\right) = \eta^{\mu}_{Y}I_{\mu}\left(t\right) }[/math]

任务，从可观测量 [math]\displaystyle{ X_{\mu}\left(t\right) }[/math]的已知时间序列中，甚至 [math]\displaystyle{ X_{\mu}\left(t\right) =\sum_{\mu}X_{\mu}\left(t\right) }[/math]中估计出来再生数^[4]（有效再生数的求解方法就是把上面的方程在不动点附近展开成线性近似，然后求出来Jacobian矩阵的本征值，当所有本征值都小于0——相当于再生数等于1——的时候，就得到了再生数。最好能够区分来自于不同接触类型[math]\displaystyle{ \beta^{\mu\nu}, \alpha^{\mu\ nu}, \gamma^{\mu}, \eta_{X,Y}^{\mu} }[/math]的有效再生数）。注意初始值[math]\displaystyle{ S_{D}\left(0\right), L_{D}\left(0\right), I_{D}\left(0\right), X_{D}\left(0\right), Y_{D}\left(0\right) }[/math]的影响，以及外生变量[math]\displaystyle{ \xi\left(t\right), s_{O}, l_{O}, i_{O} }[/math]对有效再生数的影响。

两区域简化模型：千万要减少区域间流动

考虑一个只有两个区域，区域内部均匀混合，区域之间存在流动的传染病模型。为了简化问题，我们仅仅考虑SI（假设治愈之后可以再次感染）。假设两个区域完全对称，假设流动人口不受特殊隔离。如果假设被收治或者被治愈或者被死亡之后，就不再进入接触和下一轮病情，则还需要加入[math]\displaystyle{ R }[/math]，变成SIR。

1. [math]\displaystyle{ \frac{d}{dt}S_{1}\left(t\right) = -\beta\frac{S_{1}\left(t\right)I_{1}\left(t\right)}{N_{1}\left(t\right)} +\gamma I_{1}\left(t\right) + J^{2}_{1}\frac{S_{2}\left(t\right)}{N_{2}\left(t\right)}-J^{1}_{2}\frac{S_{1}\left(t\right)}{N_{1}\left(t\right)} }[/math]
2. [math]\displaystyle{ \frac{d}{dt}I_{1}\left(t\right) = \beta\frac{S_{1}\left(t\right)I_{1}\left(t\right)}{N_{1}\left(t\right)} -\gamma I_{1}\left(t\right) + J^{2}_{1}\frac{I_{2}\left(t\right)}{N_{2}\left(t\right)}-J^{1}_{2}\frac{I_{1}\left(t\right)}{N_{1}\left(t\right)} }[/math]
3. [math]\displaystyle{ \frac{d}{dt}S_{2}\left(t\right) = -\beta\frac{S_{2}\left(t\right)I_{2}\left(t\right)}{N_{2}\left(t\right)} +\gamma I_{2}\left(t\right) + J^{1}_{2}\frac{S_{1}\left(t\right)}{N_{1}\left(t\right)}-J^{2}_{1}\frac{S_{2}\left(t\right)}{N_{2}\left(t\right)} }[/math]
4. [math]\displaystyle{ \frac{d}{dt}I_{2}\left(t\right) = \beta\frac{S_{2}\left(t\right)I_{2}\left(t\right)}{N_{2}\left(t\right)} -\gamma I_{2}\left(t\right) + J^{1}_{2}\frac{I_{1}\left(t\right)}{N_{1}\left(t\right)}-J^{2}_{1}\frac{I_{2}\left(t\right)}{N_{2}\left(t\right)} }[/math]

其中，[math]\displaystyle{ N_{\mu}\left(t\right)= S_{\mu}\left(t\right)+I_{\mu}\left(t\right) }[/math]。为了进一步简化，我们把这里的[math]\displaystyle{ N_{\mu}\left(t\right) }[/math]当做常数（实际上它们是t的函数，不过，由于将来我们只做无病平衡点附近的线性稳定性分析，这个变化是个高阶项，忽略问题不大）。

顺便，线性稳定性分析的意思是，在某个状态下，由于某个因素——例如初始条件产生了随机的小小的变化——产生了对这个状态的一个偏离，系统是否会回到这个状态，还是从这个状态离开跑到其他状态上去。在传染病里面，这个无病平衡点附近的线性稳定性分析的含义就是，看看系统对于没有被感染者的这个状态来说，一个小小的偏离——例如忽然多了很少的数量的病人——是否会造成传染病的传播。

于是，我们得到无病平衡点为[math]\displaystyle{ \left(S_{1},I_{1},S_{2},I_{2}\right)=\left(N,0,N,0\right) }[/math]。对这个状态做线性稳定性分析，也就是把原始的微分方程线性化，得到

[math]\displaystyle{ \frac{d}{dt}\left(\begin{matrix}\delta S_{1}\left(t\right)\\ \delta I_{1}\left(t\right) \\ \delta S_{2}\left(t\right) \\ \delta I_{2}\left(t\right)\end{matrix}\right) = \left[\begin{matrix}-j & -\beta+\gamma & j & 0 \\ 0 & \beta-\gamma - j & 0 & j \\ j & 0 & -j & -\beta+\gamma\\ 0 & j & 0 & \beta-\gamma-j\end{matrix}\right]\left(\begin{matrix}\delta S_{1}\left(t\right)\\ \delta I_{1}\left(t\right) \\ \delta S_{2}\left(t\right) \\ \delta I_{2}\left(t\right)\end{matrix}\right) }[/math]

其中[math]\displaystyle{ j=\frac{J}{N} }[/math]。

可以算出来，中间的矩阵的本征值为[math]\displaystyle{ \beta - 2j - \gamma, \beta - \gamma, -\sqrt{2}j, \sqrt{2}j }[/math]。我们发现，当[math]\displaystyle{ j=0 }[/math]的时候，回到已知的结果，[math]\displaystyle{ \beta - \gamma, \beta - \gamma, -0, 0 }[/math]，也就是说，当[math]\displaystyle{ \beta \gt \gamma }[/math]的时候，本征值大于零，会造成病人数量增长，也就是传染病能传播。于是，自然就得到这个模型的基本再生数为[math]\displaystyle{ R_{0}=\frac{\beta}{\gamma} }[/math]。现在，我们来看[math]\displaystyle{ j\neq 0 }[/math]的时候，我们发现其中大于零的本征值[math]\displaystyle{ \sqrt{2}j }[/math]会一直存在，也就是说，理论上，只要小小的[math]\displaystyle{ j\neq 0 }[/math]就会造成传染病的传播。当然，实际上，真的很小的时候，只需要对流动人口隔离或者专门收治，或者无差别隔离和收治，就可以抵抗这个小小的[math]\displaystyle{ j }[/math]。

但是，我们一定要看到，就这个小小的[math]\displaystyle{ j }[/math]就会把原来不会造成传染的系统[math]\displaystyle{ R_{0}=\frac{\beta}{\gamma}\lt 1 }[/math]变成一个能够传染的系统。如果愿意，我们可以在方程中引入治疗或者隔离措施的描述，看看，有和没有这个[math]\displaystyle{ j }[/math]对这些措造成的若果要保证把系统调回来到不传染状态的压力的不同。

完整的把[math]\displaystyle{ N_{\mu}\left(t\right) }[/math]看做演化变量的计算，就只能靠求微分方程的数值解了。我这里就不求了。

多个城市的反应扩散方程，考虑武汉输入和其他城市输入

只要在前面的方程上增加[math]\displaystyle{ {ij}, {\mu\nu} }[/math]角标就可以，例如[math]\displaystyle{ S^{i}_{\mu} }[/math]表示[math]\displaystyle{ i }[/math]城市的来自于[math]\displaystyle{ \mu }[/math]的易感人群数量。方程形式类似。

考虑有目标的控制措施的模型

基于多源数据的传染病计算平台及其在风险准备和风险应对中的应用

对于传染病的风险准备和风险应对来说，我们主要回答以下问题

1. 不同级别（按照传染性、死亡率、治疗成本、防控成本等分类）的传染病的发生概率大概如何；
2. 针对不同级别的传染病，各个省市或者地区的卫生应急中心，应该准备什么样的资源，做什么样的防控措施预案；
3. 在传染病实际发生阶段，如何在了解传染病的现状和预测未来，以及在此的基础上做更好的资源调配和采取更合适的防控措施。
4. 在传染病快结束的阶段，在一定风险阈值下，允许在什么地点优先放开什么防控措施。
5. 为了从数据和方法上支持以上的风险准备和风险应对，交通数据、疾病调查数据、传染病模型的建立参数估计和模拟如何调整和整合成一个计算平台，以及在现有条件下整合出来这样一个平台。

主要解决方法是

1. 运用安全学科灾害风险管理的方法，通过调查历史数据，来回答第一个关于不同级别传染病发生概率的问题
2. 以包含人口信息和经济信息的网格地理数据、交通OD矩阵（源点-终点矩阵）数据和交通流量数据、疾病传播的传染病调查数据或者设定的疾病传播模型参数，构建多源数据传染病计算平台
3. 在此计算平台上，如果是风险准备情景，则针对给定的传染病传播数据和防控应对措施，算出来传播结果，以及用于不同防控措施下的损失评估
4. 在此计算平台上，如果是传染病已经发生的风险应对情景，则结合已有的传染病传播数据，推断传染病传播参数，结合哭能的防控应对措施，算出来传播结果，以及用于不同防控措施下的损失评估

计算平台技术细节

1. 我们先来做一些设定。
   1. 区域的设定：任何一个级别的单位，例如省、市、县、地理信息系统的网格等都可以是我们计算平台上的一个区域；原则上我们的模型可以把所有的区域都等权计算，也就是如果有[math]\displaystyle{ m=1,2,\cdots, M }[/math]个这样的区域则我们的计算所用的变量关注每一个区域，例如[math]\displaystyle{ t }[/math]时刻每个区域的病人数量[math]\displaystyle{ I_{m}\left(t\right) }[/math]。但是，在传染病问题中，往往会有一个源头区域，一个所关注的目标区域。这个时候，我们就可以把所有的区域做一些组合和简化，例如分成特定源头区域[math]\displaystyle{ D }[/math]，所关注的目标区域[math]\displaystyle{ D }[/math]，以及其他区域[math]\displaystyle{ O }[/math]。这样计算模型的变量就会大大缩减。类似的，OD矩阵和交通流量数据也可以相应地合并成只有这三个区域的OD情形。当然，类似的可以把三区域模型拓展成为四区域模型等等。我们还把这样的根据具体问题做了区域合并之后的多区域模型称为[math]\displaystyle{ N }[/math]区域模型。我们把源头区域记为区域[math]\displaystyle{ 0 }[/math]。
   2. 区域间人流量的设定：根据交通OD矩阵和流量数据，以及设定的限流或者封城等防控措施，任意两区域之间的人流量[math]\displaystyle{ C^{m}_{n}\left(t\right) }[/math]是完全已知的。这里我们当做计算模型的外生变量。如果有必要，我们也可以建立非常简化的模型——例如用当前时刻之前的一个窗口的平均——来“预测”区域间人流量。
2. 模型核心计算思路：针对实际发生的传染病的风险应对情景从传染病传播的历史数据估计传染病的有效再生数参数（如果是风险准备情景，则直接给定所假设的传染病的基本再生数来当做这个传染病的一段时间的有效再生数），包含[math]\displaystyle{ R^{m}_{n}\left(t\right) }[/math]（区域 [math]\displaystyle{ m }[/math]输入到区域 [math]\displaystyle{ n }[/math]的病人的有效再生数）和 [math]\displaystyle{ R^{n}_{n}\left(t\right) }[/math]（区域 [math]\displaystyle{ n }[/math]本地传播的有效再生数）；然后，对这些有效再生数按照给定的防控措施情景或者简单预测模型——例如用当前时刻之前的一个窗口平均值来当做下一个时刻的值——来预测有效再生数；最后，利用预测出来的有效再生数，结合交通人流量来计算每个区域的病人数量。结合经济学投入产出模型做不同防控措施下的损失估计。
3. 模型：
   1. 基于有效再生数，区分了实际病人数量和观测到（被收治和统计）病人数量的模型
      1. 在[math]\displaystyle{ t }[/math]时刻从区域[math]\displaystyle{ m }[/math]输入到区域[math]\displaystyle{ n }[/math]的病人的数量，[math]\displaystyle{ I^{0,m}_{n}\left(t,0\right) }[/math]。注意，[math]\displaystyle{ I^{0,m}_{n}\left(t,0\right) \sim C^{m}_{n}\left(t\right) }[/math]和区域间人口流量密切相关
      2. 在[math]\displaystyle{ t }[/math]时刻从区域[math]\displaystyle{ m }[/math]输入到区域[math]\displaystyle{ n }[/math]的病人的数量[math]\displaystyle{ \tau }[/math]时间以后剩下的数量，[math]\displaystyle{ I^{0,m}_{n}\left(t,\tau\right)=I^{0,m}_{n}\left(t,0\right)\left(1-\sum_{\hat{\tau}=1}^{\tau}O^{0,m}_{n}\left(\hat{\tau}\right)\right) }[/math]，其中[math]\displaystyle{ O^{0,m}_{n}\left(\tau\right) }[/math]是一个从区域[math]\displaystyle{ m }[/math]输入到区域[math]\displaystyle{ n }[/math]的病人在输入之后的第[math]\displaystyle{ \tau }[/math]天被收治的概率。
      3. 在[math]\displaystyle{ t }[/math]时刻被收治的从区域[math]\displaystyle{ m }[/math]输入到区域[math]\displaystyle{ n }[/math]的病人数量，[math]\displaystyle{ X^{0,m}_{n}\left(t\right)=\sum_{\tau=1}^{\infty} X^{0,m}_{n}\left(t, \tau\right)=\sum_{\tau=1}^{\infty}I^{0,m}_{n}\left(t-\tau,0\right)O^{0,m}_{n}\left(\tau\right) }[/math]
      4. 在[math]\displaystyle{ t }[/math]时刻被从区域[math]\displaystyle{ m }[/math]输入到区域[math]\displaystyle{ n }[/math]的输入病人感染的一代感染者数量，[math]\displaystyle{ I^{1,m}_{n}\left(t,0\right)=\sum_{\tau=1}^{\infty}I^{0,m}_{n}\left(t-\tau, \tau\right)R^{0\rightarrow 1,m}_{n}\left(t\right)\omega^{0,m}_{n}\left(\tau\right) }[/math]
      5. 在[math]\displaystyle{ t }[/math]时刻被从区域[math]\displaystyle{ m }[/math]输入到区域[math]\displaystyle{ n }[/math]的病人感染的一代病人过了[math]\displaystyle{ \tau }[/math]时间以后剩下的数量，[math]\displaystyle{ I^{1,m}_{n}\left(t,\tau\right)=I^{1,m}_{n}\left(t,0\right)\left(1-\sum_{\hat{\tau}=1}^{\tau}O^{1,m}_{n}\left(\hat{\tau}\right)\right) }[/math]，其中[math]\displaystyle{ O^{1,m}_{m}\left(\tau\right) }[/math]是一个一代病人在被感染之后的第[math]\displaystyle{ \tau }[/math]天被收治的概率
      6. 在[math]\displaystyle{ t }[/math]时刻被收治的一代感染者数量[math]\displaystyle{ X^{1,m}_{n}\left(t\right)=\sum_{\tau}X^{1,m}_{n}\left(t,\tau\right) = \sum_{\tau}I^{1,m}_{n}\left(t-\tau,\tau\right)O^{1,m}_{n}\left(\tau\right) }[/math]
      7. 本地城市二代以及二代以上以上感染者数量[math]\displaystyle{ I^{2+,n}_{n}\left(t,0\right)=\sum_{\tau=1}^{\infty}\left[\sum_{m\neq n}I^{1,m}_{n}\left(t-\tau,\tau\right)+I^{2+,n}_{n}\left(t-\tau,\tau\right)\right]R^{1+\rightarrow 2+,n}_{n}\left(t\right)w^{1+,n}_{n}\left(\tau\right) }[/math]
      8. 在[math]\displaystyle{ t }[/math]时刻被本地感染的病人过了[math]\displaystyle{ \tau }[/math]时间以后剩下的数量，[math]\displaystyle{ I^{2+,n}_{n}\left(t,\tau\right)=I^{2,n}_{n}\left(t,0\right)\left(1-\sum_{\hat{\tau}=1}^{\tau}O^{2+,n}_{n}\left(\hat{\tau}\right)\right) }[/math]，其中[math]\displaystyle{ O^{2+,n}_{n}\left(\tau\right) }[/math]是一个本地二代以及二代以上病人在输入之后的第[math]\displaystyle{ \tau }[/math]天被收治的概率
      9. 在[math]\displaystyle{ t }[/math]时刻被收治的本地二代以及二代以上病人数量，[math]\displaystyle{ X^{2+,n}_{n}\left(t\right)=\sum_{\tau=1}^{\infty} X^{2+,n}_{n}\left(t, \tau\right)=\sum_{\tau=1}^{\infty}I^{2+,n}_{n}\left(t-\tau,0\right)O^{2,n}_{n}\left(\tau\right) }[/math]
      10. 为了运用这个模型，疾病调查数据需要提供[math]\displaystyle{ X^{0,m}_{n}\left(t, \tau\right), X^{1,m}_{n}\left(t, \tau\right), X^{2+,n}_{n}\left(t, \tau\right) }[/math]，每一个零代病人到达本区域之后到被收治的时间，每一个一代病人从被传染（和零代接触）之后到被收治的时间，每一个二代病人从被传染（和一代病人接触）之后到被收治的时间。有了这些数据，就能估计出来零代病人传染力函数[math]\displaystyle{ \omega^{0,m}_{n}\left(\tau\right) }[/math]、一代以及一代以上病人传染力函数[math]\displaystyle{ w^{1+,n}_{n}\left(\tau\right) }[/math]、零代一代二代及以上病人的收治概率函数[math]\displaystyle{ O^{0,m}_{n}\left(\tau\right), O^{1,m}_{n}\left(\tau\right), O^{2+,n}_{n}\left(\tau\right) }[/math]，然后进一步用来估计出来从零代感染一代的有效再生数[math]\displaystyle{ R^{0\rightarrow 1,m}_{n}\left(t\right) }[/math]、从一代以及以上感染二代以及以上的有效再生数[math]\displaystyle{ R^{1+\rightarrow 2+,n}_{n}\left(t\right) }[/math]。更进一步，基于现有[math]\displaystyle{ R^{0\rightarrow 1,m}_{n}\left(t\right), R^{1+\rightarrow 2+,n}_{n}\left(t\right) }[/math]的简单预测，以及外生变量输入人口[math]\displaystyle{ C^{m}_{n}\left(t\right) }[/math]以及来源城市[math]\displaystyle{ m }[/math]内的病人密度和病程分布，就可以用于预测病人数量。
      11. 模型局限性以及应对：这个模型假设收治概率函数[math]\displaystyle{ O\left(\tau\right) }[/math]仅仅是病程[math]\displaystyle{ \tau }[/math]的函数而不依赖于被收治时间点[math]\displaystyle{ t }[/math]。实际上，它应该是[math]\displaystyle{ t,\tau }[/math]两者的函数。可是，真的都变成两者的函数，整个模型的参数推断部分就做不出来了（请数学高手再看看是否可以做出来），于是，我们不得不暂时接受这个假设。另外，确实一个社会整体收治能力的变化是一个远远比病程引起的收治概率来的慢。那么，可以怎么来超越这个局限性呢？可以取一个时间窗口的数据来估计[math]\displaystyle{ O\left(\tau\right) }[/math]。然后，在这个时间窗口后面一小段时间，采用这个统计出来[math]\displaystyle{ O\left(\tau\right) }[/math]。实际计算过程中，也可以不断地做时间窗口的平移，直到发现一个比较大的改变点，然后必要的时候采用后面算出来的[math]\displaystyle{ O\left(\tau\right) }[/math]。
   2. 区分了实际病人数量和观测到（被收治和统计）病人数量的微分方程模型（这部分的参数估计和估计完之后用于预测，还有待解决；微分方程模型和上面的实证模型的结合和比较，也是一个要解决的问题）
   3. 用于估算防控措施效果的多主体模型（结合交通数据、经济-人口地理网格数据、出行方式和频率的改变、接触方式和频率的改变）
   4. 三种模型（传染病时间序列建模、传染病微分方程模型、传染病多主体模型）的综合和比较：针对同一个假设的或者已经发生的传染病的数据，把不同模型给出的结果做一个对比；如果发现具有相容性和互补性，必要的时候模型之间可以相互结合——一个模型的某些数据进入另一个模型来相互补充。

计算平台创新之处

1. 风险管理分析方法、交通OD数据和流量数据、疾病传播调查数据、疾病传播模型的参数估计和模拟、经济学投入产出分析等多源数据-模型-分析方法的结合
2. 在传染病模型的参数拟合和模拟中，区分观测到病人数量和实际病人数量。在之前的建模和分析中，研究者大多数假设实际病人数量是直接可测的，或者说在实证研究参数拟合中倒过来把观测到的病人数量直接当做实际病人数量。基于这个新的模型，我们已经和原作者合作一起改进了传染病有效再生数计算软件EpiEstim。
3. 在传染病模型的参数拟合和模拟中，区分输入病人的有效再生数[math]\displaystyle{ R^{1,m}_{n}\left(t\right) }[/math]和本地病人的有效再生数[math]\displaystyle{ R^{2,n}_{n}\left(t\right) }[/math]，尽管对疾病传播调查提出了新的数据需求，但是，能够做到更加准确的计算。基于这个新的区分，我们已经和原作者合作一起改进了传染病有效再生数计算软件EpiEstim。
4. 系统科学的系联思想的独特作用：不能就一个区域或者一个系统（例如仅仅考虑传染病本身）来考虑，必须考虑联系，联系有传染、交通、资源、经济生产和消费等多种，需要综合考虑。

参考文献