一个由规则定义的系统：使用元胞自动机模拟并预测传染病传播

抽象

新型冠状病毒(COVID-19)的爆发对全球经济和人们的生活都造成了严重的影响．在本篇文章中我们使用元胞自动机模型和遗传算法模拟并预测病毒的传播．该模型具有简单易懂，易实现，参数少，参数易学习的特点，同时，可结合遗传算法使用．我们只需少量的外部数据，因为我们的模型可以高效地从数据中挖掘出有效信息．一个元胞自动机的系统是由规则定义的，我们只需声明出该系统的规则，系统就会自动地模拟出这种规则下事物的演化方式和相互作用的迭代过程，使得用元胞自动机对离散系统进行仿真和模拟非常简单.

关键词

新型冠状病毒 传染病 元胞自动机 遗传算法

基本情况介绍

新型冠状病毒（以下简称新冠）于2019年年尾已经有零星的确诊案例被报道，于2020年农历新年左右爆发，该病毒具有极强的传染性，并且潜伏期较长，目前为止没有特效药可用于治疗，亦无有效的疫苗，戴上医用外科口罩或防护性较好的N95(KN95，或等同的标准）口罩能大幅减少被传染的风险，而中国政府快速响应实施全民居家隔离措施，得以快速控制疫情，可见隔离还是主要的应对新冠的有效方法，而足够的检测和饱和的医疗支援供应对于疫情的控制来说也非常重要.

元胞自动机模型和人员流动表示

我们首先将人群划分为三类：健康者记做h(Health的首字母)，无症状感染者也叫潜伏者记做c(Conceal的首字母)，确诊感染者记做i(Infected的首字母)．然后我们用xm表示带上口罩的人(m取自Mask的首字母)，不管是健康者还是潜伏者还是感染者，即xm可能是hm也可能是cm也可能是im，例如，其中hm表示带上口罩的健康者，我们再用xq表示被隔离的人(q取自Quarantine首字母)，用xc表示正在接受治疗的人(c取自Cure首字母).

我们用a代表空气(Air)，也就是空间，可让人自由移动.

人们所处的空间可看成是一条收尾相连的环，这条环上有一个个连接着的格子，人们可在格子间移动，如下图，我们定义了一个4个格子的空间，一个健康者，要表示一个健康者4天在一直走动，可以这么表示，为简单起见我们假设只有一个方向的移动：

day1: a h a a
day2: a a h a
day3: a a a h 
day4: h a a a
day5: a h a a
...

在2维的元胞自动机中，一条规则可以这么表示

_ _ _
  _

或者记做

(_, _, _) -> _

例如

| h a _ | _ h a |
|   h   |   a   |

或者等价地，可以记做

(h, a, _) -> h
(_, h, a) -> a

表示：对所有健康患者，如果他前面是空气，那么他移动，并且，它原来的位置变成空气，而它移动到的那个位置要从空气变为他．其实就是，如果匹配到了 h a _ 这样的模式，那么中间的a要替换为h，而如果匹配到了 _ h a 这样的模式，那么中间的h要替换为a．你可以拿这个规则代入到day1-day5那个例子去验证．其中的下划线表示，不管那个格子里写的是什么符号都匹配.

隔离患者是不可以随便移动的，因此有

(_, hq, _) -> hq

为了减少人员流动，可以定义一个自发隔离的概率

(h, a, _) -> h, (p=p1)
(h, a, _) -> hq, (p=p2)
p1 + p2 = 1

其中表示，有一定的概率，一个健康人移过来之后立刻被隔离，这样它就不能再移动了，被隔离的概率是p1．显然p1越大，人员流动率越小．人员流动率可以用p2表示, p2=1-p1．

元胞自动机的运作的实质无非就是模式的匹配和替换，用白话来说，就是「如果我怎么样，如果我的左邻居是怎么样，如果我的右邻居是怎么样，那么我接下来怎么样」，例如，

h, a, _ -> h

表示：如果我的左邻居是健康者，那么我被替换成健康者，而不管我的右邻居是什么.

_, h, a -> a

表示：如果我的右邻居是空气，而不管我的左邻居是什么，那么我变为空气，这两条规则放在一起

h, a, _ -> h
_, h, a -> a

就实现了一个健康者h在空气中的不断移动．

类似这样地，整个系统所有的运作方式都可以用一系列的规则来表示：

规则的构成的格式其实就是模式 -> 替换值，例如，我们来看这条规则：

(c|i)m?, h, _ -> c

它的模式是

(c|i)m?, h, _

右邻居是什么，不影响，但是它要求自己是h，然后自己的左邻居或者是c，或者是i，而且戴不戴口罩都无所谓，以下任何一个模式，都是这条模式的子模式：

c, h, _
i, h, _
cm, h, _
im, h, _

也就是说

pattern1: (c|i)m?, h, _

这条模式实际上就是

pattern2: c, h, _
pattern3: i, h, _
pattern4: cm, h, _
pattern5: im, h, _

这四条模式写在一起，我们给这些模式标上了记号pattern几，pattern1被满足，当且仅当pattern2到pattern5这四条pattern的任何一条被满足．其他的模式，如果你觉得写法奇怪，它们都是类似的道理.

例如说，如果元胞自动机的状态是

a a c h h

那么由于

c h h

匹配了

c, h, _

这个模式，所以，中间的h会被替换成c，元胞自动机的状态就下一次迭代中演变成

a a c c h

类似地，假如有

a i h h h

这样的模式，由于

i h h

能够匹配

i h _

这个模式，所以根据规则，它中间的h要被替换成c，该元胞自动机的下一次的状态就是

a i c h h

接着根据规则，c又会传染右边的h，这个h变成c后又相继再传染给最右边的h，这个元胞自动机接下来的状态就是

a i c c h
a i c c c

这样，传染病的传播就在元胞自动机中得到了表示.

总的思路

元胞自动机本身其实不会给出预测，一个定义良好的元胞自动机描述的仅仅是怎么从第i轮迭代的状态状态i演化到第i+1轮迭代的状态状态i+1:

CA: 状态i --> 状态i+1

如果我们想要用元胞自动机来做预测，我们得知道初始状态，也就是说有初始状态和元胞自动机，才能做预测，然后其实就是元胞自动机的一次次演化就会演化到接近现今结果的状态．我们说元胞自动机模型，就是把元胞自动机看成是一个数学模型，类似线性回归模型，线性回归模型有参数要估计，而同样地，我们建出来的这个元胞自动机模型，它也有参数要估计，这个参数是什么呢？你已经猜到了，这个参数就是初始状态.

模型的参数估计就是要把这个模型的参数找出来，使得模型的准确率最高，参数在哪里找呢？在搜索空间里面找．例如，我们在估计一元线性回归模型

$$ y = \alpha + \beta x + \mu, \quad \mu \approx N(0, \sigma^2) $$

的回归方程

$$ y = \hat{\alpha} + \hat{\beta}x + e $$

的时候我们实际上是在找$\alpha$的估计$\hat{\alpha}$和$\beta$的估计$\hat{\beta}$使得残差

$$ e = y - \hat{y} = y - (\hat{\alpha} + \hat{\beta}x) $$

的平方和

$$ \sum e^2 $$

变得尽可能的小，也就是当所有其他的$\hat{\alpha}$和$\hat{\beta}$的取值都不能使$\sum e$变得更小的时候，$\hat{\alpha}$和$\hat{\beta}$就算是找到了.

抽象地来看，这里有一个损失函数$\sum e^2$它受参数$(\hat{\alpha}, \hat{\beta})$的影响，模型参数估计的过程，实际上就是找到使损失函数最小的参数的过程．也就是说，要找到这个模型的参数，我们首先得为这个模型定义一个损失函数，在一元线性回归模型中，我们定义的损失函数是关于$\hat{\alpha}$和$\hat{\beta}$的残差平方和函数．

元胞自动机模型的参数，我们已经知道了，是初始状态，这是我们要去确定的值，有了这个才能用元胞自动机去推演和预测．我们如果把时间周期定为天，我们把元胞自动机系统里的每一次迭代与现实生活中的每一天的时间流逝对应起来，那么，假设每一天的疫情数据是

$$ y_1, y_2, y_3, y_4, \cdots $$

而设元胞自动机的初始状态也是$y_1$，那么元胞自动机就可以演化出结果

$$ y_1, \hat{y_2}, \hat{y_3}, \hat{y_4}, \cdots $$

这时，损失函数可以合理地定义为

$$ \sum (y_i - \hat{y_i})^{2} $$

其中的$\hat{y_i}$实际上就表示元胞自动机模型第$i-1$次迭代得到的值，直观来说就是，如果这个元胞自动机模型是好的，是比较准确的，那它至少能比较接近地解释过去发生了的事情，在这个基础之上，我们才能开始尝试拿它去预测未来．比如说，第1天，第2天，第3天，第4天的累积确诊人数和累积出院人数是这样的

$$ (100, 1), (120, 2), (200, 3), (310, 4) $$

那么，至少应该找到一种初始状态使得元胞自动机输出的是比较接近的值，例如

$$ (99, 0), (132, 4), (198, 5), (330, 4) $$

元胞自动机的输出结果也是状态，怎么得到累积确诊人数和累积治愈人数呢？记得我们用i表示确诊的元胞，我们统计匹配模式i的元胞的数量就可以了.

在Mathematica中实现它

Mathematica自带了很强大的元胞自动机函数

CellularAutomaton[规则，初始状态，步数]

基于它，我们可以很方便地实现我们要实现的模型.

疾病传播的简单演示

上图是一个用Mathematica的元胞自动机来模拟疾病传播的演示过程，首先我们定义符号i表示感染者，h表示健康者，a表示空气，当h旁边有i时，h会变成i，当h和i隔着a时，h不会变成i，这是最简单的模型，它对应的规则集是

i, h, _ -> i
_, h, i -> i
_, i, _ -> i;
h|a, h, h|a -> h;
_, a, _ -> a

它对应的Mathematica代码是

Clear[f]

f[{i, h, _}] := i;
f[{_, h, i}] := i;
f[{_, i, _}] := i;
f[{h | a, h, h | a}] := h;
f[{_, a, _}] := a;

MatrixForm[CellularAutomaton[
    {f[#] &, {}, 1},
    {h, h, i, h, h, h, a, a, a, h, h, h, a},
    5
]]

其中的

{h, h, i, h, h, h, a, a, a, h, h, h, a}

其实是初始状态，我们可以看到，矩阵的一行表示一次迭代的状态，左边的人群有一个感染者，而很快，感染就扩散到整个人群，右边的人群和左边的人群隔着a，因此没有被传染.

下面，我们来模拟进行了一定程度的交通管控后，人们的出行情况如何，首先，我们要定义一个由一系列规则组成的规则集：

_, h, a -> a
h, a, _ -> maybe[h, hq, 0.3]
_, a, _ -> a
_, hq, _ -> hq
_, x_, _ -> x
_, hq, _ -> maybe[h, hq, 0.7]

其中，这两个maybe就减慢了人们的流动速度，因为maybe[a,b,p]实际上表示「取a或者b，而取到a的概率为p」．maybe[a,b,p]的函数定义如下

maybe[x_, y_, p_] := (
    RandomVariate[
        BinomialDistribution[1, p]
    ] /. {1 -> x, 0 -> y}
);

而模拟这个交通管控的元胞自动机的Mathematica代码如下：

Clear[f]

f[{_, h, a}] := a;
f[{h, a, _}] := maybe[h, hq, 0.3];
f[{_, a, _}] := a;
f[{_, hq, _}] := hq;
f[{_, x_, hq}] := x;
f[{_, hq, _}] := maybe[h, hq, 0.7];
f[{_, h, h}] := h;

ArrayPlot[
    CellularAutomaton[
        { f[#] &, {}, 1 },
        {a, a, a, a, h, a, a, a, a, a, a, a, a},
        10
   ] /. hq -> h, 
   ColorRules -> {h -> White, a -> Black}
]

在Mathematica中模拟交通管控

可以看到交通管控之后h（白色的）的行走速度是下降了，因为有那么几次迭代h一步都没有．即便是有多个人员，系统仍能正常运行：

system = CellularAutomaton[
    { f[#] &, {}, 1 },
    {a, h, a, a, h, a, a, a, h, a, a, a, a},
    60
];

ArrayPlot[
    system /. hq -> h, 
    ColorRules -> {h -> White, a -> Black}
]

模拟交通系统仍能稳定运行

当然我们并不仅仅满足于用元胞自动机模拟一个人员能够受控流动的交通系统，虽然说病毒的传播需要人员流动因而交通系统的模拟是整个系统及其重要的组成，但是，我们最初的目的是要模拟和预测传染病的传播，因此，接下来还要在实现中引入感染者.

colorRules = {
  h -> White,
  a -> Black,
  i -> Red
};

symbolSubstitutes = {
   hq -> h,
   iq -> i
};

Clear[f];

f[{h, a, _}] := maybe[h, hq, 0.3];
f[{h, i, _}] := maybe[i, iq, 0.3];
f[{_, iq, _}] := maybe[i, iq, 0.7];
f[{_, hq, _}] := maybe[h, hq, 0.7];

f[{_, h, a}] := a;
f[{_, a, _}] := a;
f[{_, hq, _}] := hq;
f[{_, x_, hq}] := x;
f[{_, h, h}] := h;
f[{_, i, _}] := i;
f[{i | (iq), h, _}] := i;
f[{_, h, i | (iq)}] := i;
f[{i, a, _}] := i;
f[{_, i, a}] := a;

initialStates = {
    a, a, h, h, a, h, a, a, h, h, a, a, h, a, i, a, a
};

system = CellularAutomaton[
    {f[#] &, {}, 1}, 
    initialStates,
    120
];

ArrayPlot[
    system /. symbolSubstitutes, 
    ColorRules -> colorRules
]

元胞自动机模型引入了感染者

这个模型最初只有1位感染者，演进了120步之后，大部分的人都已经被感染，正如图中所示的，白色的是健康者，而红色的是感染．进一步地，我们可以统计累积感染者的数量变化：

infectedCounts = #[i] & /@ (
    Counts[#] & /@ (
        system /. symbolSubstitutes
    )
);

ListPlot[infectedCounts, Filling -> Axis]

输出结果如下图所示

感染者总人数走势

可以看到，由于实施了交通管控，并且空余的空间较多，所以感染人数增加得比较慢，病毒传播得比较慢.

parabola = Fit[infectedCounts, {1, x, x^2}, x];

Show[
    ListPlot[infectedCounts, Filling -> Axis], 
    Plot[parabola, {x, 0, 120}]
]

呈二次曲线增长

趋势线接近一条二次多项式的曲线．现在有一个问题是，一个2维的元胞自动机究竟能不能模拟感染人数指数增长的情况呢？我们的答案并不是太确定，考虑到新冠具有潜伏期，当人员流动速度加快，当无症状感染者流动至各个地方，当多个地点同时爆发，即使是在2维的元胞自动机中，也可能会有亚指数增长的现象出．

即使，这种2维的元胞自动机模型模拟不出指数增长的情形，我们也仍然会花时间来研究它，这是因为，任何一个模型它都有自身的局限性，天下没有免费的午餐，不会存在这么一种万能的通用的模型能够用来解决一切问．况且，就我国的疫情曲线来看，确诊人数也未必应被认为是指数增长的.

模型的基本假设

如上所述，任何模型都有其自身的局限性和适用范围，我们建出来的「元胞自动机」模型也不例外， 在这样一个一维的元胞自动机中，感染数量的增长速度不会发展得太快，并且，每一个具有传染性的元胞也只能把病毒传染给它旁边的两个元胞，由此，我们可以总结出这个模型所要求的两条假定：

1. 病毒的传播速度具有自限性，或者存在这么一种传播「阻力」，当病毒传播速度变快时，「阻力」也变大，如同下落的物体的速度不会一直增加，而是会因为阻力最终达到一个恒定值一．就新冠病毒而言，当疫情发展到一定程度，民众的卫生意识会觉醒，政府也会开始行动，区别仅仅是行动的时间是早或晚而已.

2. 病毒是在固定范围传播的，假设A携带病毒，B离A最近，C离B最近，而C不是离A最近的点，这样，A会首先传播给B，然后B才传播给C，而不是A同时传播给B和C，并且，实际上，A会同时传播给B和C，当且仅当，B和C都是离A最近的点.

感染人数可能会在前期快速增加，但是到了一定阶段，增加的占比最终会减小，体现为增速放缓， 而累积感染人数的曲线也会变得平缓，这一点，和实际是比较吻合的，因此，这两条假定一定程度上来说是合理的.

模型的进一步演进

接下来我们要引入一个随机的初始状态生成器，因为，想要看看更大规模的系统的演化过程是什么样的，以及，即使确定了确诊人数，还是不能确定元胞自动机的初始状态，因为，有若干种元胞自动机的状态都能算出同样的确诊人数，换言之，对应于一个确诊人数的元胞自动机的状态是不唯一的，因此元胞自动机的初始状态可以看成我们想要拟合的参数，下面我们从初始状态的搜索空间中随机找到一个点，即产生一个随机的初始状态：

randomInitialStates[
   healthies_Integer, 
   infecteds_Integer, 
   conceals_Integer,
   airs_Integer
] := (
    RandomSample[
        Join[
            Table[h, healthies],
            Table[i, infecteds],
            Table[c, conceals],
            Table[a, airs]
        ]
    ]
);

Mathematica这种声明式的语言，使用起来真的非常简便：

生成随机的初始状态

下面为使我们的模型更进一步，我们修改规则，使健康人h被感染后不立刻成为确诊者i，而是先成为携带者c，c再依概率演化为确诊者i.

colorRules = {
  h -> White,
  a -> Black,
  i -> Red
};

symbolSubstitutes = {
   hq -> h,
   iq -> i,
   c -> h,
   cq -> h
};

Clear[f];

param1 = 0.3;
param2 = 0.3; 
param3 = 0.3;
param4 = 0.3;
param5 = 0.3;
param6 = 0.3;
param7 = 0.3;
param8 = 0.3;
param9 = 0.3;

f[{h, a, _}] := maybe[h, hq, param1];
f[{_, hq, _}] := maybe[h, hq, param2];

f[{i, a, _}] := maybe[i, iq, param3];
f[{_, iq, _}] := maybe[i, iq, param4];

f[{c, a, _}] := maybe[c, maybe[cq, iq, param8], param5];
f[{_, cq, _}] := maybe[c, maybe[cq, iq, param9], param6];

f[{_, h, a}] := a;
f[{_, i, a}] := a;
f[{_, c, a}] := a;
f[{_, c, Except[a]}] := maybe[i, c, param7];

cannotMove = (a | iq | hq | cq);
f[{cannotMove, a, _}] := a;
f[{_, h, h}] := h;
f[{_, i, Except[a]}] := i;

infectious = ( i | iq | c | cq);

f[{
    infectious,
    h,
    _
}] := c;

f[{
    _,
    h,
    infectious
}] := c;

f[{
    Except[(i | iq | c | cq)], 
    h, 
    h | hq
}] := h;

f[{
    Except[(i | iq | c | cq)],
    hq,
    h | hq
}] := hq;

initialHealthies = 7;
initialInfecteds = 2;
initialConceals = 6;
initialAirs = 10;
initialStates = randomInitialStates[
    initialHealthies,
    initialInfecteds,
    initialConceals,
    initialAirs
];

stepsToEvolve = 80;

system = CellularAutomaton[
    {f[#] &, {}, 1}, 
    initialStates,
    stepsToEvolve
];

ArrayPlot[
    system /. symbolSubstitutes, 
    ColorRules -> colorRules
]

基于上述规则，当系统的初始状态为：

{
    a, i, a, h, a, c, a, c, h, h,
    c, c, i, a, h, h, a, a, a, c,
    c, a, h, h, a
}

的时候，系统演化80步的演化过程如下图所示：

系统演化80步

以下代码画出感染人数的走势图：

healthyCounts[system_] := (
    #[h] & /@ (
        Counts[#] & /@ (system /. symbolSubstitutes)
    )
);

infectedCounts[system_] := (
    #[i] & /@ (
        Counts[#] & /@ (system /. symbolSubstitutes)
    )
);

ListPlot[infectedCounts[system], Filling -> Axis]

输出结果如下图：

感染人数走势图

可以看到，在这个系统中，前期的感染人数的增加是非常快的，这是因为有很多分散各地的无症状感染者，而到后期感染人数增长速度则放缓，直到最后所有人都被感染.

接下来我们需要引入治愈措施，使得i或者iq可以以一定的概率变为h，为此我们需要更改规则：

colorRules = {
  h -> White,
  a -> Black,
  i -> Red
};

symbolSubstitutes = {
   hq -> h,
   iq -> i,
   c -> h,
   cq -> h
};

Clear[f];

param1 = 0.3;
param2 = 0.3; 
param3 = 0.3;
param4 = 0.3;
param5 = 0.3;
param6 = 0.3;
param7 = 0.3;
param8 = 0.3;
param9 = 0.3;
param10 = 0.3;
param11 = 0.3;
param12 = 0.3;

f[{h, a, _}] := maybe[h, hq, param1];
f[{_, hq, _}] := maybe[h, hq, param2];

f[{i, a, _}] := maybe[i, maybe[iq, h, param11], param3];
f[{_, iq, _}] := maybe[i, maybe[iq, h, param12], param4];

f[{c, a, _}] := maybe[c, maybe[cq, iq, param8], param5];
f[{_, cq, _}] := maybe[c, maybe[cq, iq, param9], param6];

f[{_, h, a}] := a;
f[{_, i, a}] := a;
f[{_, c, a}] := a;
f[{_, c, Except[a]}] := maybe[i, c, param7];

cannotMove = (a | iq | hq | cq);
f[{cannotMove, a, _}] := a;
f[{_, h, h}] := h;
f[{_, i, Except[a]}] := maybe[h, i, param10];

infectious = ( i | iq | c | cq);

f[{
    infectious,
    h,
    _
}] := c;

f[{
    _,
    h,
    infectious
}] := c;

f[{
    Except[(i | iq | c | cq)], 
    h, 
    h | hq
}] := h;

f[{
    Except[(i | iq | c | cq)],
    hq,
    h | hq
}] := hq;

initialHealthies = 7;
initialInfecteds = 2;
initialConceals = 6;
initialAirs = 10;
initialStates = randomInitialStates[
    initialHealthies,
    initialInfecteds,
    initialConceals,
    initialAirs
];

stepsToEvolve = 80;

system = CellularAutomaton[
    {f[#] &, {}, 1}, 
    initialStates,
    stepsToEvolve
];

引入治愈率之后我们再来看系统的演化过程：

ArrayPlot[system /. symbolSubstitutes, ColorRules -> colorRules]

会发现这时感染人数最终会得到控制，感染比例也是被控制着：

感染比例被控制着

感染人数的走势如图：

感染人数走势

下边我们引入全新的符号系统来应对日益复杂的状态表示和演算需求：

原子：

标头和复合标头：

Clear[
    p, f, proceed, dispatch
];

Clear[
    movable, unMovable, infectious, nonInfectious,
    veryInfectious, lessInfectious,
    veryInfecttable, lessInfecttable, infecttable, nonInfecttable
];

p[x_Integer] := 0.3;

movable = ( a | h | i | c | _m | _t );
unMovable = ( _s | _st | _sm );

lessInfectious = (
    sm[i] |
    sm[c] |
    m[i] |
    m[c]
);

veryInfectious = (
    i |
    c |
    s[i] |
    s[c]
);

infectious = (
    lessInfectious |
    veryInfectious
);

nonInfecttable = (
    _t |
    _st |
    infectious |
    a
);

lessInfecttable = (
    m[h] |
    sm[h]
);

veryInfecttable = (
    h |
    s[h]
);

infecttable = (
    lessInfecttable |
    veryInfecttable
);

nonInfectious = (
    h |
    m[h] |
    sm[h] |
    _t |
    _st |
    a |
    s[h]
);

proceed[x:a] := (
    a
);

proceed[x:h] := (
    RandomChoice[
        { p[1], p[2], p[3], p[4] } -> {
            h, m[h], s[h], sm[h]
        }
    ]
);

proceed[x:i] := (
    RandomChoice[
        { p[5], p[6], p[7], p[8], p[9], p[10] } -> {
            s[i], t[i], st[i], a, sm[i], m[i]
        }
    ]
);

proceed[x:c] := (
    RandomChoice[
        { p[11], p[12], p[13], p[14], p[15], p[16], p[17] } -> {
            i, m[i], sm[i], c, m[c], s[c], sm[c]
        }
    ]
);

proceed[x_s] := (
    RandomChoice[
        { p[26], p[27] } -> {
            x, proceed[x[[1]]]
        }
    ]
);

proceed[x_m] := (
    RandomChoice[
        { p[28], p[29] } -> {
            x, proceed[x[[1]]]
        }
    ]
);

proceed[x_sm] := (
    RandomChoice[
        { p[30], p[31] } -> {
            x, proceed[x[[1]]]
        }
    ]
);

proceed[x_t] := (
    RandomChoice[
        { p[18], p[19], p[20], p[21], p[25] } -> {
            m[h], sm[h], x, a, st[x[[1]]]
        }
    ]
);

proceed[x_st] := (
    RandomChoice[
        { p[22], p[22], p[23], p[24], p[40] } -> {
            m[h], sm[h], x, a, t[x[[1]]]
        }
    ]
);

infectFromSource[_, x:nonInfecttable] := x;

infectFromSource[source:nonInfectious, y_] := y;

infectFromSource[source:veryInfectious, y:veryInfecttable] := (
    RandomChoice[
        { p[32], p[33] } -> {
            y, c
        }
    ]
);

infectFromSource[source:veryInfectious, y:lessInfecttable] := (
    RandomChoice[
        { p[34], p[35] } -> {
            y, m[c]
        }
    ]
);

infectFromSource[source:lessInfectious, y:veryInfecttable] := (
    RandomChoice[
        { p[36], p[37] } -> {
            y, c
        }
    ]
);

infectFromSource[source:lessInfectious, y:lessInfecttable] := (
    RandomChoice[
        { p[38], p[39] } -> {
            y, m[c]
        }
    ]
);

infectFromLeftSide[leftSideX_, y_] := (
    infectFromSource[leftSideX, y]
);

infectFromRightSide[x_, rightSideY_] := (
    infectFromSource[rightSideY, x]
);

infectFromDoubleSource[x_, y:nonInfecttable, z_] := y;

infectFromDoubleSource[x_, y:infecttable, z_] := (
    RandomChoice[
        { p[41], p[42] } -> {
            y, c
        }
    ]
);

infectFromBothSide[x_, y_, z_] := (
    infectFromDoubleSource[x, z, y]
);

proceed[x:infectious, y_, z:nonInfectious] := (
    infectFromLeftSide[x, y]
);

proceed[x:nonInfectious, y_, z:infectious] := (
    infectFromRightSide[y, z]
);

proceed[x:infectious, y_, z:infectious] := (
    infectFromBothSide[x, y, z]
);

proceed[x:nonInfectious, y_, z:nonInfectious] := (
    proceed[y]
);

dispatch[x:movable, a, _] := (
    proceed[x]
);

dispatch[_, x:movable, a] := (
    a
);

dispatch[x:unMovable, a, _] := (
    a
);

dispatch[_, x:unMovable, a] := (
    proceed[x]
);

dispatch[x_, y_, z:Except[a]] := (
    proceed[x, y, z]
);

dispatch[x_, y:Except[a], z_] := (
    proceed[x, y, z]
);

f[{x_, y_, z_}] := dispatch[x, y, z];

initialHealthies = 5;
initialInfecteds = 2;
initialConceals = 3;
initialAirs = 16;
initialStates = randomInitialStates[
    initialHealthies, 
    initialInfecteds, 
    initialConceals, 
    initialAirs
];

stepsToEvolve = 2;

system = CellularAutomaton[
    {f[#] &, {}, 1}, 
    initialStates,
    stepsToEvolve
];

displaySystem[system_] := (
    symbolSubstitute = {
        s[h] -> h,
        sm[h] -> h,
        _t -> i,
        _st -> i,
        c -> h,
        s[c] -> h,
        sm[c] -> h,
        s[i] -> i,
        sm[i] -> i,
        m[i] -> i,
        m[h] -> h,
        m[c] -> h
    };

    colorRules = {
        h -> White,
        i -> Red,
        a -> Black
    };

    ArrayPlot[
        system /. symbolSubstitute,
        ColorRules -> colorRules
    ]
);

在某种参数下情况最终得到控制

为了简便，我们设定了默认参数

p[_Integer] = 0.3;

在这种条件下，即使考虑到无症状感染者，考虑到戴口罩的人，考虑到收治情况，情况最终还是得到了控制（如图．那么现在我们有了一个基本完整的元胞自动机，给定了初期的感染人数，不仅有元胞的布局需要确定，还要确定那四十多个参数，确定这些参数的规则是，使得对历史的与预测值与实际值相差最小.

数据的收集

我们用的数据呢，是从互联网上搜索得到的数据，具体是参考这一个网页的数据，这份数据是知乎网友爬取相关的疫情数据网站，然后整理得到．第一列表示日期，格式是月份-日数，第二列表示累积确诊数，第三列表示疑似病例，第四列表示累积治愈，第五列表示累积死亡人．我们现在用的这份数据可以从[这里下载](/covid-19-data.csv．这一份数据描述的是国内的1月12号到3月18号的感染情况.

使用遗传算法训练和拟合模型

观察上边我们写出来的元胞自动机代码可以知道，有p[1]一直到p[42]这42个参数待确定（我们使用的是缺省值），并且，即使知道累积有多少人感染，初始状态也没有确定（我们使用的是随机生成的一个初始状态），观察疫情数据，我们大概可以确定的是，有两个阶段，前一个阶段，无症状携带者很多，而确诊的仅仅是少部分，因此增速越来越快，后面的阶段，大部分患者都已经被收治（感谢国家的「应收尽收」政策），无症状携带者也越来越少，因此呢，累积确诊数的新增也就越来越小，这条疫情曲线，初期到中期非常陡，到中后期和后期，就越来越平缓．我们提到这两个阶段，是由于这个元胞自动机模拟的系统，它的疫情发展速度是近似是线性的，所以，前面那个阶段，我们要单独训练一个元胞自动机，后面那个阶段，又要另外训练一个元胞自动机，不同的阶段，对应（或者说适用于）不同的模型.

中国大陆1月12日到3月18日的新冠确诊人数变化图

看上面这幅图，在转折点之前累积感染人数是加速上升的，而在转折点之后，累积感染人数的增加速度（新增确诊人数）逐渐减小，这个转折点大概是在2月13号左．这样我们就可以用1月12日到2月13日的数据训练一个模型，然后再用2月14日到3月18日的数据训练另外一个模．实际上，出于实用考虑，我们只关注用2月14日到3月18的数据训练出来的模型就足够了.

在遗传算法的框架中有几个重要的概念：

1. 染色体
2. 表现型
3. 表现型得分

遗传算法主要是通过「变异」（类似基因突变）随机改变染色体中的某个碱基位点引入新的基因型，表现型实际上则是我们要求的解的备选，在这个案例中，表现型就是我们要确定的元胞自动机的参数，而表现型得分也就是适应度，是这个表现型，也就是这个备选解它的优劣评估，例如一个好的表现型（一列好的参数）总能使模型至少较为准确地解释历史事件，这时该表现型的得分就高一些，而一个不好的表现型使模型给出的理论预测值与历史观测偏差很大，使模型甚至连历史事件都解释不了，更加不可能用来预测未来事件了，因此，这样的表现型的表现型得分肯定是很低，遗传算法会使得分低的表现型对应的染色体被抛弃，使得分高的表现型对应的染色体留存下．类似自然的演化规律.

总结起来，遗传算法不断重复的无非是这几个步骤：

1. 筛选适应度高的个体留存下来
2. 留存下来的个体参与繁衍产生下一代
3. 随机地引入基因突变
4. 跳到步骤1重复直至找到满意的表现型

对于我们的这个元胞自动机模型来说，表现型是那42个权重参数，外加一列元胞初始状态，而在给定了表现型之后，得分是由模型的预测值和历史数据的观测值的偏差平方和确定的，即，得分是由数据确定的，代入什么样的数据，就学习出什么样的模．难的是确定这样一种从表现型到染色体的转换规则：假设我的表现型是x，我怎样得到与之对应的染色体y？假设染色体是y，怎样反算出对应的表现型x？

观察我们的元胞符号，只有a h c i这四个是最基本的，而s[_] m[_] t[_] sm[_] st[_]实际上都是系统演化过程中产生的，因此，我们不需要把带标头的符号像s[_] m[_] .．st[_]这些纳入到表现型的考虑中，因为元胞自动机的初始状态可以没有它们（它们在演化过程中出现），要注意到，数据显示到最后有八万多人感染，但是我们的元胞自动机模型的初始状态并不会有那么多元胞，这实际上是涉及到了一个「颗粒度」的问题，设想要模拟一个化学反应，可以在分子团的尺度去模拟，模拟分子团和分子团之间的相互作用，还可以在分子的尺度模拟，模拟分子和分子之间的相互作用，甚至还可以在原子的尺度模拟，模拟原子和原子之间的相互作用，但是，到了原子尺度，原子的数量是非常多的，所以计算机要维护和计算的基本单元的数量也非常多，计算负担会变得难以承受，我们的元胞自动机模型也是这个道理，我们不会在人的尺度去模拟，我们最终实际上不会让一个元胞对应到一个人，而是，让一个元胞对应到一个基本的可相互作用的单元，例如，h可以对应一个易受感染的单元，这样的一个单元可以代表好几百人，m[h]可以对应到一个加强防护单元，每个单元对应几百个人，这样子，每一轮迭代要处理的元胞数量大大减小，计算变得真正可行了，到最后我们只需比较元胞自动机模型得到的疫情曲线和历史数据的疫情曲线的相似度就可以了，如果这两条曲线很贴近，我们就说这个模型是具有一定的解释能力的，因而也能参考它的预测输．而对于那42个参数，简单起见，就不用再抽象了，42个参数可以不用再精简了.

表现型到染色体的转换过程图示

要由变现型得到染色体（如图），首先分别编码参数部分和元胞自动机初始状态，对参数部分进行编码得到一列二进制数，再对元胞自动机初始状态进行编码，也得到一列二进制数，将这两列二进制数首位相接，得到一列拼接的二进制数，这一个拼接得到的二进制数，就是「染色体」，就可以用于遗传算法中的交叉和变．而把染色体转化会表现型也是同样的过程反过来，我们要记住那个区分参数编码和元胞初始状态编码的位点，例如参数编码出了n位二进制数，而元胞初始状态编码出了m位二进制数，那么，拼接二进制数中的第1到第m位就对应参数的编码，而第m+1到第m+n位对应元胞初始状态的编码，有了参数编码可以转换回一个个参数，有了元胞初始状态编码也可以转换回元胞初始状态，编解码过程大概就是这样子的.

一个参数转化为二进制数是这样的：

encodeParameter[parameter_] := (
    PadLeft[IntegerDigits[parameter-1, 2], 4]
);

对于一个参数列表，它被编码为二进制数的过程是这样的

encodeParameters[parameters_] := (
    Flatten @ ( 
        encodeParameter[#] & /@ parameters
    )
);

一个元胞转换成二进制编码的过程是这样的

encodeCell[cell_] := (
    cellToNum = {a -> 0, h -> 1, i -> 2, c -> 3};
    num = cell /. cellToNum;
    PadLeft[IntegerDigits[num, 2], 2]
);

一列元胞转换成二进制编码的过程是这样的

encodeCells[cells_] := (
    [email protected](encodeCell[#] & /@ cells)
);

给定初始状态和初始元胞状态，要得到二进制编码，转换过程是这样的：

encodeStates[parameters_, cells_] := (
    Join[
        encodeParameters[parameters],
        encodeCells[cells]
    ]
);

得到了一个（一组）参数的编码，要转换回原来的一个（一组）参数的过程是这样的：

decodeToParameter[code_] := FromDigits[code, 2] + 1;

decodeToParameters[code_] := (
    decodeToParameter /@ ArrayReshape[
        code, { Length[code]/4, 4 }
    ]
);

类似地，得到了一个（一组）元胞的编码，要转换会原来的一个（一组）元胞，过程是这样的：

decodeToCell[code_] := FromDigits[code, 2] /. {
    0 -> a,
    1 -> h,
    2 -> i,
    3 -> c
};

decodeToCells[code_] := decodeToCell /@ ArrayReshape[
    code, { Length[code]/2, 2 }
];

知道了一个编码蕴含了多少个参数和多少个元胞，我们也可以由这个编码得到原来的参数和元胞组：

decodeToPhenotype[code_, nParameters_Integer, nCells_Integer] := (
    parameters = decodeToParameters[
        Take[code, { 1, nParameters*4 }]
    ];

    cells = decodeToCells[Take[
        code, 
        { 4*nParameters + 1, 4*nParameters + 2*nCells }
    ]];

    Association[{"parameters" -> parameters, "cells" -> cells}]
);

事实上，可以很容易地验证这些转换规则是否工作正常：

验证转换规则是否正确

可以看到编码功能和解码功能都能正常工．

遗传算法用到的初始种群一般而言都是随机的，或者是上一次遗传算法运行得到的，我们这里呢，也是先写一个函数，这个函数随机地产生nParameters个参数，nCells个元胞组成的初始状态：

randomPhenotype[nParameters_, nCells_] := (
    parameters = Table[
        RandomChoice[
            Table[i, {i, 1, 16}]
        ],
        nParameters
    ];

    cells = Table[
        RandomChoice[
            {a, h, c, i}
        ],
        nCells
    ];

    Association[{
        "parameters" -> parameters,
        "cells" -> cells
    }]
);

以及生产随机初始种群的函数

randomChromosome[nParameters_, nCells_] := (
    phenotype = randomPhenotype[nParameters, nCells];
    encodeStates[phenotype["parameters"], phenotype["cells"]]
);

再把元胞自动机也做成函数，做成一个「演化算子」的形式：

cellularAutomatonMaker[parameters_] := (

    p = parameters;

    movable = ( a | h | i | c | _m | _t );
    unMovable = ( _s | _st | _sm );

    lessInfectious = (
        sm[i] |
        sm[c] |
        m[i] |
        m[c]
    );

    veryInfectious = (
        i |
        c |
        s[i] |
        s[c]
    );

    infectious = (
        lessInfectious |
        veryInfectious
    );

    nonInfecttable = (
        _t |
        _st |
        infectious |
        a
    );

    lessInfecttable = (
        m[h] |
        sm[h]
    );

    veryInfecttable = (
        h |
        s[h]
    );

    infecttable = (
        lessInfecttable |
        veryInfecttable
    );

    nonInfectious = (
        h |
        m[h] |
        sm[h] |
        _t |
        _st |
        a |
        s[h]
    );

    proceed[x:a] := (
        a
    );

    proceed[x:h] := (
        RandomChoice[
            { p[1], p[2], p[3], p[4] } -> {
                h, m[h], s[h], sm[h]
            }
        ]
    );

    proceed[x:i] := (
        RandomChoice[
            { p[5], p[6], p[7], p[8], p[9], p[10] } -> {
                s[i], t[i], st[i], a, sm[i], m[i]
            }
        ]
    );

    proceed[x:c] := (
        RandomChoice[
            { p[11], p[12], p[13], p[14], p[15], p[16], p[17] } -> {
                i, m[i], sm[i], c, m[c], s[c], sm[c]
            }
        ]
    );

    proceed[x_s] := (
        RandomChoice[
            { p[26], p[27] } -> {
                x, proceed[x[[1]]]
            }
        ]
    );

    proceed[x_m] := (
        RandomChoice[
            { p[28], p[29] } -> {
                x, proceed[x[[1]]]
            }
        ]
    );

    proceed[x_sm] := (
        RandomChoice[
            { p[30], p[31] } -> {
                x, proceed[x[[1]]]
            }
        ]
    );

    proceed[x_t] := (
        RandomChoice[
            { p[18], p[19], p[20], p[21], p[25] } -> {
                m[h], sm[h], x, a, st[x[[1]]]
            }
        ]
    );

    proceed[x_st] := (
        RandomChoice[
            { p[22], p[22], p[23], p[24], p[40] } -> {
                m[h], sm[h], x, a, t[x[[1]]]
            }
        ]
    );

    infectFromSource[_, x:nonInfecttable] := x;

    infectFromSource[source:nonInfectious, y_] := y;

    infectFromSource[source:veryInfectious, y:veryInfecttable] := (
        RandomChoice[
            { p[32], p[33] } -> {
                y, c
            }
        ]
    );

    infectFromSource[source:veryInfectious, y:lessInfecttable] := (
        RandomChoice[
            { p[34], p[35] } -> {
                y, m[c]
            }
        ]
    );

    infectFromSource[source:lessInfectious, y:veryInfecttable] := (
        RandomChoice[
            { p[36], p[37] } -> {
                y, c
            }
        ]
    );

    infectFromSource[source:lessInfectious, y:lessInfecttable] := (
        RandomChoice[
            { p[38], p[39] } -> {
                y, m[c]
            }
        ]
    );

    infectFromLeftSide[leftSideX_, y_] := (
        infectFromSource[leftSideX, y]
    );

    infectFromRightSide[x_, rightSideY_] := (
        infectFromSource[rightSideY, x]
    );

    infectFromDoubleSource[x_, y:nonInfecttable, z_] := y;

    infectFromDoubleSource[x_, y:infecttable, z_] := (
        RandomChoice[
            { p[41], p[42] } -> {
                y, c
            }
        ]
    );

    infectFromBothSide[x_, y_, z_] := (
        infectFromDoubleSource[x, z, y]
    );

    proceed[x:infectious, y_, z:nonInfectious] := (
        infectFromLeftSide[x, y]
    );

    proceed[x:nonInfectious, y_, z:infectious] := (
        infectFromRightSide[y, z]
    );

    proceed[x:infectious, y_, z:infectious] := (
        infectFromBothSide[x, y, z]
    );

    proceed[x:nonInfectious, y_, z:nonInfectious] := (
        proceed[y]
    );

    dispatch[x:movable, a, _] := (
        proceed[x]
    );

    dispatch[_, x:movable, a] := (
        a
    );

    dispatch[x:unMovable, a, _] := (
        a
    );

    dispatch[_, x:unMovable, a] := (
        proceed[x]
    );

    dispatch[x_, y_, z:Except[a]] := (
        proceed[x, y, z]
    );

    dispatch[x_, y:Except[a], z_] := (
        proceed[x, y, z]
    );

    f[{x_, y_, z_}] := dispatch[x, y, z];

    evolver = CellularAutomaton[ {f[#] &, {}, 1} ];

    evolver
);

那么既然编码器和解码器都正常工作了，我们来试一下，首先我们生成一组随机的参数，并且用这组参数定义一个元胞自动机：

phenotype = randomPhenotype[42, 100];

parameters = AssociationThread[
    Table[i, {i, 1, 42}] -> phenotype["parameters"]
];
cells = phenotype["cells"];

ca = cellularAutomatonMaker[parameters];

然后我们再用制造出来的这个元胞自动机来演化一个初始状态：

system = NestList[ca, cells, 40];
displaySystem[system]

结果如图所示：

用元胞自动机演化一个随机生成的初始状态

可以看到疫情最后还是被消灭．

下一步呢，我们就是要找到这么一组参数，使得元胞自动机在这种配置下，疫情的发展情况能和现实中的类似，怎么找呢？当然是用遗传算法了，遗传算法某种程度上来说呢，是通用的，遗传算法你可以看做是定义在染色体域上的一种计算，这也就意味着遗传算法的代码是可以复用的.

rData = Drop[data, 1];
unCureds = rData[[All, 2]] - rData[[All, 4]] - rData[[All, 5]];

infectedCounts[system_] := (

    Clear[i];

    symbolSubstitutes = {
        s[h] -> h, 
        sm[h] -> h, 
        _t -> i, 
        _st -> i, 
        c -> h,
        s[c] -> h, 
        sm[c] -> h, 
        s[i] -> i, 
        sm[i] -> i, 
        m[i] -> i, 
        m[h] -> h, 
        m[c] -> h
    };

    #[i] & /@ (
        Counts[#] & /@ (system /. symbolSubstitutes)
    ) /. _Missing -> 0
);

computeScoresForSystem[system_] := (
    systemOutputs = infectedCounts[system];
    N[
        Correlation[
            systemOutputs,
            Take[unCureds, -Length[systemOutputs]]
        ],
        12
    ]
);

generatePopulation[nPopulation_Integer] := (
    nParameters = 42;
    nCells = 400;
    phenotypes = Table[randomPhenotype[nParameters, nCells], nPopulation];

    encodeStates[
        #["parameters"], #["cells"] 
    ]& /@ phenotypes
);

rawNumberListToParameters[numList_] := (
    AssociationThread[Table[i, {i, 1, 42}] -> numList]
)

computeScoresForChromosome[chromosome_] := (

    (* must have exact same value as what generatePopulation does *)
    nParameters = 42;
    nCells = 400;
    nStepsToEvolve = 20;

    phenotype = decodeToPhenotype[chromosome, nParameters, nCells];
    caParameters = phenotype["parameters"];
    caCells = phenotype["cells"];
    caMachine = cellularAutomatonMaker[
        caParameters // rawNumberListToParameters
    ];
    system = NestList[caMachine, caCells, nStepsToEvolve];

    scores = computeScoresForSystem[system];
    scores
);

randomBitsFlips[chromosome_] := (
    flipRate = 0.09;
    flipOrNot = RandomVariate[
        BinomialDistribution[1, flipRate], 
        Length[chromosome]
    ];

    newChromosome = Table[
        If[
            flipOrNot[[i]] == 1, 
            RandomInteger[{0, 1}], 
            chromosome[[i]]
        ],
        {i, 1, Length[chromosome]}
    ]
)

populationVary[populationx_] := (
    varyRate = 0.11;
    varyOrNot = RandomVariate[
        BinomialDistribution[1, varyRate], 
        Length[populationx]
    ];
    varied = Table[
        If[
            varyOrNot[[i]] == 1, 
            randomBitsFlips[populationx[[i]]], 
            populationx[[i]]
        ],
        {i, 1, Length[populationx]}
    ]
)

sampleCross[x_, y_] := (

    midPoint[chromosome_ /; EvenQ[Length[chromosome]]] := (
        Length[chromosome] / 2
    );

    midPoint[chromosome_ /; OddQ[Length[chromosome]]] := (
        (Length[chromosome]-1) / 2
    );

    {
        Join[
            Take[x, {1, midPoint[x]}], 
            Take[y, {midPoint[x]+1, Length[y]}]
        ],

        Join[
            Take[y, {1, midPoint[x]}], 
            Take[x, {midPoint[x]+1, Length[y]}]
        ]
    }
);

populationCross[populationx_] := (
    populationMale = Take[
        populationx, 
        {1, Length[populationx] - 1}
    ];
    
    populationFemale = Take[
        populationx, 
        {2, Length[populationx]}
    ];

    crossed = Table[
        sampleCross[populationMale[[i]], populationFemale[[i]]],
        {i, 1, Length[populationMale]}
    ];

    ArrayReshape[
        Flatten[crossed], 
        {
            Dimensions[crossed][[1]]*Dimensions[crossed][[2]], 
            Length[populationx[[1]]]
        }
    ]
)

dropUnfits[populationScores_] := (
    dropRate = 0.64;
    populationSize = Length[populationScores];
    selectScores = TakeSmallest[
        populationScores, 
        Round[populationSize*(1 - dropRate)]
    ]
)

filterPopulation[population_] := (
    dropRate = 0.64;
    sortedPopulation = 
    ReverseSortBy[population, computeScoresForChromosome];
    Take[sortedPopulation, 
    Round[(1 - dropRate)*Length[sortedPopulation]]]
);

gaEvolver[population_] := (
    Clear[selectedPopulations, variedPopulations, nextGeneration];
    selectedPopulations = filterPopulation[population];
    variedPopulations = randomBitsFlips /@ selectedPopulations;
    nextGeneration = populationCross[variedPopulations];
    Join[variedPopulations, nextGeneration]
);

定义好遗传算法演算子之后，就可以开始演算了，首先我们生成一个规模为100的初始种群，然后让该种群演化10代，看每一代的最高得分者的得分表现：

initialPopulation = generatePopulation[100];
generations = NestList[gaEvolver, initialPopulation, 10];
highestScores = Table[
    computeScoresForChromosome[generation[[1]]],
    { generation, generations }
];
ListLinePlot[highestScores]

大多数点都位于0.95以上

这里的得分，就是每一代里面表现最好的那个元胞自动机的预测结果和历史数据的相关系数，得分先是趋于增长，后趋于稳定，我们可以选表现最优的的那一个元胞自动机：

lastGeneration = Last[generations];
bestChromosome = First[lastGeneration];
bestPhenotype = decodeToPhenotype[bestChromosome, 42, 400];
bCaParameters = bestPhenotype["parameters"];
bCaCells = bestPhenotype["cells"];
bCaMachine = cellularAutomatonMaker[
    bCaParameters // rawNumberListToParameters
];
bNStepsToEvolve = 20;
bSystem = NestList[bCaMachine, bCaCells, bNStepsToEvolve];
bSystemOutputs = infectedCounts[bSystem];
ListPlot[
    {
        bSystemOutputs*300, 
        Take[unCureds, -bNStepsToEvolve]
    }, 
    PlotLabels -> {"predict", "historical"}
]

预测值对比历史值

元胞自动机中模拟的系统的感染人数走势，和历史数据的感染人数走势的曲线相似，这正是我们想要的效果，这个元胞自动机模型是用历史数据的后20个值来拟合的，如果我们想要向前预测几个值，即，历史数据最终是到3月18号，我们想看3月19号，3月20号，或者3月21号的值，我们只需让元胞自动机模型向前演化23步，因为演化20步对应3月18号的数据，21步对应19号的，22步对应20号的，而23步就对应21号的，来试试吧：

predictSystem = NestList[bCaMachine, bCaCells, 23];
predictSystemOuts = infectedCounts[predictSystem];

ListPlot[
    {
        predictSystemOuts*300, 
        Take[unCureds, -bNStepsToEvolve]
    }, 
    PlotLabels -> {"predict", "historical"}
]

changeRate1 = predictSystemOuts[[21]]/predictSystemOuts[[20]];
changeRate2 = predictSystemOuts[[22]]/predictSystemOuts[[20]];
changeRate3 = predictSystemOuts[[23]]/predictSystemOuts[[20]];

predictUncureds1 = Round[N[Last[unCureds]*changeRate1]]
predictUncureds2 = Round[N[Last[unCureds]*changeRate2]]
predictUncureds3 = Round[N[Last[unCureds]*changeRate3]]

画出来的图是这样的：

蓝色的点是预测值的走势

3月19号到21号的未治愈感染人数预测值如下：

8183
8664
9146

写了一个星期的任务，终于完成辣.

总结

对于新新冠状病毒疫情的预测，我们考虑了两个假设：

1. 病毒的传播具有自限性，当疫情严重到一定程度之后，居民和政府会开始行动，居民卫生意识增强，医疗资源强化供应，公共出行减少，病毒传播速度减缓.

2. 病毒的传播具有区域性，A是B最近的点，B是C最近的点，C不是A最近的点，A携带病毒，那么A首先将病毒传播给B，然后B再将病毒传播给C，A会将病毒同时传播给B和C，当且仅当，B和C都是A最近的点.

在这两个假设下，我们用一个1维的元胞自动机模型来模拟病毒的传播机制，一开始我们模拟人的行走，然后我们模拟受管控的交通，然后我们模拟病毒的无阻力传播，然后我们加入更复杂的参数，模拟具有医疗救治和防护的情形，我们将现实中的单元簇抽象为元胞自动机里的一个元胞，也就是说一个元胞可以对应到多个人，然后考虑元胞自动机演化过程中感染人数的变化趋势，如果这个变化趋势和现实中得到的历史数据的变化相似，我们就采用这个模型，就可以用这个模型做预．我们用了遗传算法来寻找模型的参数，使得模型在这组参数下和历史数据拟合得尽可能的好，并且遗传算法也成功的使拟合效果越来越．元胞自动机模型不能模拟所有的病毒传播，尤其是当现实中的病毒传播太快时，元胞自动机模型跟不上，所以，它适用于模拟疫情发展后期的情况.

参考文献

[1] CellularAutomaton—Wolfram Language Documentation

[2] Coronavirus disease 2019 - Wikipedia