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

背景

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

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]

1. 在[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]
2. 在[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]天被收治的概率
3. 在[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]

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

如何从实际数据得到[math]\displaystyle{ O_{D}\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]，先求比例在平均和先加总（平均）再求比例。

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

新冠肺炎相关计算

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

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

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{ D_{t} }[/math]：武汉输入病人数量的存量。假设[math]\displaystyle{ t\leq 0 }[/math]的时候封城，之后不再有新来的武汉病人，都是封城之前积累的病人，总数记为[math]\displaystyle{ D_{0} }[/math]。之后每天随着被发现和死亡减少。假设武汉输入病人的传染性可以用每天的传染性参数[math]\displaystyle{ k_{D} }[/math]表示。暂时我们假设这个参数不依赖于时间。如果有必要，可以考虑这个参数的时间序列。
2. [math]\displaystyle{ d_{t} }[/math]：每天确诊的武汉输入病人的数量，流量。假设[math]\displaystyle{ d_{t} }[/math]（[math]\displaystyle{ t\gt 0 }[/math]）数据完全已知。
3. [math]\displaystyle{ o_{t} }[/math]：每天来自于其他城市（非本地，非武汉）输入病人的数量。假设[math]\displaystyle{ o_{t} }[/math]（[math]\displaystyle{ t\gt 0 }[/math]）数据完全已知。可以通过其他城市来的人口数量乘以平均感染率（这个要通过从数据中用某种方式估算出来，或者假设一个参数[math]\displaystyle{ \mathcal{o} }[/math]加入模型）来估计。这个[math]\displaystyle{ o_{t} }[/math]可以是一个矢量，分开考虑各个城市。
4. [math]\displaystyle{ O_{t} }[/math]：来自于其他城市（非本地，非武汉）输入病人数量的存量。假设来自于其他城市的输入病人的传染性可以用每天的传染性参数[math]\displaystyle{ k_{O} }[/math]表示。原则上我们还需要一个[math]\displaystyle{ O_{0} }[/math]来表示初始的来自于其他城市的输入病人数量。不过，可以把这个初始人数算到后面的[math]\displaystyle{ M_{0} }[/math]里面。
5. [math]\displaystyle{ m_{t} }[/math]：每天本地人——通过接触武汉病人、其他城市流动过来的病人、本地已经被感染的病人——被感染得病的人的数量，流量。假设[math]\displaystyle{ m_{t} }[/math]不能直接可观测。
6. [math]\displaystyle{ M_{t} }[/math]：本地被传染得病的未收治未死亡的人的存量。假设[math]\displaystyle{ M_{0} }[/math]是数据获取起点之前就累积的病人数量。整个[math]\displaystyle{ M_{t} }[/math]，包含[math]\displaystyle{ M_{0} }[/math]不可直接观测。假设这些本地被感染有没有被发现和收治的病人的传染性可以用再生数[math]\displaystyle{ k_{3} }[/math]表示。为了简化，可以假设[math]\displaystyle{ k_{C}=k_{O} }[/math]。注意[math]\displaystyle{ k_{D} }[/math]由于管控可能会不太一样。
7. [math]\displaystyle{ x_{t} }[/math]：每天所有被传染得病的人当中被发现和收治的数量，流量。平均来说，一个本地病人每天被医疗机构发现得病的概率为 [math]\displaystyle{ p }[/math]。
8. [math]\displaystyle{ X_{t} }[/math]：累积所有被传染得病的人当中被发现和收治的数量，存量。假设被医疗机构发现和收治之后的病人的传染性可以用传染性指数[math]\displaystyle{ k_{X} }[/math]表示。为了简化，可以假设[math]\displaystyle{ k_{X}=0 }[/math]。
9. [math]\displaystyle{ y_{t} }[/math]：每天没有被发现和收治但是死亡（或者完全自我康复，也可以把死亡和自我康复两种情况分开）的人数，流量。平均来说，一个在这座城市的病人每天死亡并且没有被医疗机构发现的概率为[math]\displaystyle{ q }[/math]。假设病人死亡之后不再会造成传染。如果需要考虑死亡之后病毒的扩散，也可以引入死亡之后每天的传染性参数[math]\displaystyle{ k_{Y} }[/math]。为了简化，可以假设[math]\displaystyle{ k_{Y}=0 }[/math]。
10. [math]\displaystyle{ Y_{t} }[/math]：累积没有被发现和收治但是死亡的人数，存量。
11. 由于潜伏期也可能具有传染性，把潜伏期和传染期当做一体，合起来称为传染期。传染期在这里被当做从被感染到被发现和收治之间的一段时间。假设这段时间病人的死亡概率和传播能力和病人相同。这是一个简化假设。如果需要更细节，可以考虑一个依赖于病程时间点的传染性和死亡率。假设目标城市输入性病人传染期为[math]\displaystyle{ \tau_{D} }[/math]，其他城市输入性病人传染期为[math]\displaystyle{ \tau_{O} }[/math]，本地病人传染期为[math]\displaystyle{ \tau_{C} }[/math]。配合从个案数据获得的平均一个病人所传染的其他病人的数量[math]\displaystyle{ R_{D}, R_{C}, R_{O} }[/math]，我们可以得到每一天这三种病人的传染能力可以表达为，[math]\displaystyle{ k_{D,C,O}=\frac{\ln{\left(R_{D,C,O}\right)}}{\tau_{D,C,O}}-1 }[/math]。同时，可以通过这个[math]\displaystyle{ \tau_{D,C,O} }[/math]估计出来，特定城市输入性病人、本地病人、其他城市输入性病人的[math]\displaystyle{ p }[/math]和[math]\displaystyle{ q }[/math]。例如，[math]\displaystyle{ \tau_{C}=\frac{1}{p} }[/math]。

迭代方程

1. [math]\displaystyle{ d_{t} }[/math]，完全已知，外生变量
2. [math]\displaystyle{ D_{t}=D_{0}\left(1-q\right)^{t}-\sum_{1\leq \tau \leq t}d_{\tau} }[/math]
3. [math]\displaystyle{ o_{t} }[/math]，完全已知，外生变量
4. [math]\displaystyle{ O_{t}=\sum_{1\leq \tau \leq t}o_{\tau}\left(1-p-q\right)^{\tau} }[/math]
5. [math]\displaystyle{ m_{t}=D_{t-1}k_{D}+O_{t-1}k_{O}+M_{t-1}k_{C}+X_{t-1}k_{X}+y_{t-1}k_{Y} }[/math]，或者[math]\displaystyle{ Y_{t-1}k_{Y} }[/math]，或者[math]\displaystyle{ Y_{t-1}k_{Y}+Y_{t-2}k_{Y} }[/math]，表示这些没有被收治的死亡病人不能马上得到处理
6. [math]\displaystyle{ M_{t}=M_{t-1}\left(1-p-q\right)+m_{t}+M_{0}\left(1-p-q\right)^{t-1} }[/math]，当期不被收治
7. [math]\displaystyle{ x_{t}=d_{t}+m_{t-1}p+o_{t-1}p+M_{0}\left(1-p-q\right)^{t-1}p }[/math]，可直接观测。这里有一个简化假设：所有病人都在被感染那天以[math]\displaystyle{ p }[/math]的概率被发现和收治。
8. [math]\displaystyle{ X_{t}=\sum_{1\leq \tau \leq t}x_{\tau} }[/math]，可直接观测
9. [math]\displaystyle{ y_{t}=m_{t-1}q+O_{t-1}q+M_{0}\left(1-p-q\right)^{t-1}q }[/math]。这里有一个简化假设：所有病人都在被感染那天以[math]\displaystyle{ q }[/math]的概率死亡。
10. [math]\displaystyle{ Y_{t}=\sum_{1\leq \tau \leq t}y_{\tau} }[/math]

需要把这个模型的[math]\displaystyle{ x_{t} }[/math]（[math]\displaystyle{ X_{t} }[/math]）或者[math]\displaystyle{ x_{t} }[/math]（[math]\displaystyle{ X_{t} }[/math]）的每一项解开，然后，和实际数据对比拟合出来待定参数。

简化情况1

忽略没有得到治疗的死亡率[math]\displaystyle{ q }[/math]，也就没有了[math]\displaystyle{ Y }[/math]。同时忽略其他城市的流入，忽略医院内感染，则

1. [math]\displaystyle{ d_{t} }[/math]，完全已知，外生变量
2. [math]\displaystyle{ D_{t}=D_{0}-\sum_{1\leq \tau \leq t}d_{\tau} }[/math]
3. [math]\displaystyle{ m_{t}=D_{t-1}k_{D}+M_{t-1}k_{C} }[/math]
4. [math]\displaystyle{ M_{t}=M_{t-1}\left(1-p\right)+m_{t}+M_{0}\left(1-p\right)^{t-1} }[/math]，当期不被收治
5. [math]\displaystyle{ x_{t}=d_{t}+m_{t-1}p+M_{0}\left(1-p\right)^{t-1}p }[/math]，可直接观测
6. [math]\displaystyle{ X_{t}=\sum_{1\leq \tau \leq t}x_{\tau} }[/math]，可直接观测

只要给定[math]\displaystyle{ M_{0}, D_{0} }[/math]结合上面从实际数据统计出来的参数，其[math]\displaystyle{ X_{t} }[/math]就完全确定了。

接着忽略特定城市的流入（相当于直接计算特定城市的情况）

1. [math]\displaystyle{ m_{t}=M_{t-1}k_{C} }[/math]
2. [math]\displaystyle{ M_{t}=M_{t-1}\left(1-p\right)+m_{t}+M_{0}\left(1-p\right)^{t-1} }[/math]，当期不被收治
3. [math]\displaystyle{ x_{t}=m_{t-1}p+M_{0}\left(1-p\right)^{t-1}p }[/math]，可直接观测
4. [math]\displaystyle{ X_{t}=\sum_{1\leq \tau \leq t}x_{\tau} }[/math]，可直接观测

只要给定[math]\displaystyle{ M_{0} }[/math]结合上面从实际数据统计出来的参数，其[math]\displaystyle{ X_{t} }[/math]就完全确定了。

简化情况2

当[math]\displaystyle{ p=1 }[/math]（也就是自然[math]\displaystyle{ q=0 }[/math]了），也就是发现和收治能力足够，得病必然被发现和收治，的时候，方程简化为

1. [math]\displaystyle{ d_{t} }[/math]，完全已知，外生变量
2. [math]\displaystyle{ D_{t}=D_{0}-\sum_{1\leq \tau \leq t}d_{\tau} }[/math]
3. [math]\displaystyle{ o_{t} }[/math]，完全已知，外生变量
4. [math]\displaystyle{ O_{t}=\sum_{1\leq \tau \leq t}o_{\tau} }[/math]
5. [math]\displaystyle{ m_{t}=D_{t-1}k_{1}+O_{t-1}k_{2}+M_{t-1}k_{3} }[/math]
6. [math]\displaystyle{ M_{t}=m_{t} }[/math]
7. [math]\displaystyle{ x_{t}=d_{t}+m_{t-1}+o_{t-1} }[/math]，可直接观测。
8. [math]\displaystyle{ X_{t}=\sum_{1\leq \tau \leq t}x_{\tau} }[/math]，可直接观测

进一步，保留输入病人[math]\displaystyle{ d_{t} }[/math]，不考虑输入病人[math]\displaystyle{ o_{t} }[/math]，有

1. [math]\displaystyle{ d_{t} }[/math]，完全已知，外生变量
2. [math]\displaystyle{ D_{t}=D_{0}-\sum_{1\leq \tau \leq t}d_{\tau} }[/math]
3. [math]\displaystyle{ m_{t}=D_{t-1}k_{1}+M_{t-1}k_{3} }[/math]
4. [math]\displaystyle{ M_{t}=m_{t} }[/math]
5. [math]\displaystyle{ x_{t}=d_{t}+m_{t-1} }[/math]，可直接观测。
6. [math]\displaystyle{ X_{t}=\sum_{1\leq \tau \leq t}x_{\tau} }[/math]，可直接观测

注意：如果要用这个模型，则从[math]\displaystyle{ x_{t}=d_{t}+m_{t-1} }[/math]和[math]\displaystyle{ d_{t} }[/math]中得到[math]\displaystyle{ m_{t-1} }[/math]，然后，从[math]\displaystyle{ m_{t}=D_{t-1}k_{1}+m_{t-1}k_{3} }[/math]中可以得到[math]\displaystyle{ k_{1},k_{3} }[/math]。当然，更一般的情况，需要考虑到有限的[math]\displaystyle{ p }[/math]。

更进一步，不考虑输入病人[math]\displaystyle{ d_{t}, o_{t} }[/math]，有

1. [math]\displaystyle{ m_{t}=M_{t-1}k_{3} }[/math]
2. [math]\displaystyle{ M_{t}=m_{t} }[/math]
3. [math]\displaystyle{ x_{t}=m_{t-1} }[/math]，可直接观测
4. [math]\displaystyle{ X_{t}=\sum_{1\leq \tau \leq t}x_{\tau} }[/math]，可直接观测

于是，[math]\displaystyle{ x_{t}=M_{0}k_{3}^{t-1} }[/math]，[math]\displaystyle{ X_{t}=M_{0}\frac{1-k_{3}^{t}}{1-k_{3}} }[/math]。回到最简单的再生数模型。

这个不考虑输入病人[math]\displaystyle{ d_{t}, o_{t} }[/math]的简化，如果假设所有的病人都继续传播，而不是仅仅那些没有收治的在传播，则

1. [math]\displaystyle{ m_{t}=M_{t-1}k_{3} }[/math]
2. [math]\displaystyle{ M_{t}=M_{t-1}+m_{t}=\left(1+k_{3}\right)M_{t-1} }[/math]
3. [math]\displaystyle{ x_{t}=m_{t-1} }[/math]，可直接观测
4. [math]\displaystyle{ X_{t}=\sum_{1\leq \tau \leq t}x_{\tau} }[/math]，可直接观测

于是，[math]\displaystyle{ x_{t}=M_{0}\left(1+k_{3}\right)^{t-1}k_{3} }[/math]，[math]\displaystyle{ X_{t}=M_{0}\left(1+k_{3}\right)^{t}=M_{t} }[/math]。当然，实际上总要考虑收治（可能有限概率收治）和死亡的。

考虑潜伏期

假设潜伏期为感染后的L天，L天之内每一天的被发现概率为[math]\displaystyle{ p^{1},p^{2}, ..., p^{L} }[/math]，死亡率为[math]\displaystyle{ q^{1},q^{2}, ..., q^{L} }[/math]，自己好了的概率为[math]\displaystyle{ r^{1},r^{2}, ..., r^{L} }[/math]，其在潜伏期的感染能力用再生数[math]\displaystyle{ KL^{1},KL^{2}, ..., KL^{L} }[/math]表示。一个病人在L天之后的状态可能是被发现和收治了，死亡了，自己好了，没有被发现也没有死也没有自己好——发病了——能够继续传染。最后一种可能的几率，平均来说，是[math]\displaystyle{ 1-p^{1}-q^{1}-r^{1},1-p^{2}-q^{2}-r^{2}, ..., 1-p^{L}-q^{L}-r^{L} }[/math]。把这些概率加到迭代方程中就行：这时候需要把病人数量分组，成为处于潜伏时间为[math]\displaystyle{ l=1,2,...L }[/math]的病人的数量。

需要的时候在把这个方程写下来。

多城市的情形

1. [math]\displaystyle{ d\left(t\right)_{i} }[/math]，完全已知，外生变量
2. [math]\displaystyle{ D\left(t\right)_{i}=D\left(0\right)_{i}\left(1-q_{i}\right)^{t}-\sum_{1\leq \tau \leq t}d\left(\tau\right)_{i} }[/math]
3. [math]\displaystyle{ o\left(t\right)^{j}_{i} }[/math]，完全已知，外生变量
4. [math]\displaystyle{ O\left(t\right)_{i}=\sum_{j,1\leq \tau \leq t}o\left(\tau\right)^{j}_{i}\left(1-p_{i}-q_{i}\right) }[/math]，其他城市不包含武汉，不包含所关心的城市自己
5. [math]\displaystyle{ m\left(t\right)_{i}=D\left(t-1\right)_{i}k_{1,i}+O\left(t-1\right)_{i}k_{2,i}+M\left(t-1\right)_{i}k_{3,i}+X\left(t-1\right)_{i}k_{4,i}+y\left(t-1\right)_{i}k_{5,i} }[/math]
6. [math]\displaystyle{ M\left(t\right)_{i}=\sum_{1\leq \tau \leq t-1}m\left(\tau\right)_{i}\left(1-p_{i}-q_{i}\right)+m\left(t\right)_{i}+M\left(0\right)_{i}\left(1-p_{i}-q_{i}\right)^{t-1} }[/math]，当期不被收治，相当于假设潜伏期一个单位时间
7. [math]\displaystyle{ x\left(t\right)_{i}=d\left(t\right)_{i}+m\left(t-1\right)_{i}p_{i}+o\left(t-1\right)_{i}p_{i}+M\left(0\right)_{i}\left(1-p_{i}-q_{i}\right)^{t-1}p }[/math]，可直接观测
8. [math]\displaystyle{ X\left(t\right)_{i}=\sum_{1\leq \tau \leq t}x\left(t\right)_{i} }[/math]，可直接观测
9. [math]\displaystyle{ y\left(t\right)_{i}=m\left(t-1\right)_{i}q_{i}+O\left(t-1\right)_{i}q_{i}+M\left(0\right)\left(1-p_{i}-q_{i}\right)^{t-1}q_{i} }[/math]。这里有一个简化假设：所有病人都在被感染那天以[math]\displaystyle{ q }[/math]的概率死亡。
10. [math]\displaystyle{ Y\left(t\right)_{i}=\sum_{1\leq \tau \leq t}y\left(t\right)_{i} }[/math]

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

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

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]的易感人群数量。方程形式类似。

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

参考文献