本研究获得北京航空航天大学第三十一届“冯如杯”学生科技竞赛数理类二等奖。

本研究获得北京航空航天大学第三十一届“冯如杯”学生科技竞赛数理类二等奖。

利用 EM 算法估计混合三参数 logistic 模型的参数

Birnbaum¹ 通过对传统 logistic 模型添加参数，创造了三参数 logistic 模型，其形式如下：

其导数为：

混合三参数 logistic 模型是两个 logistic 模型的线性叠加，且每个模型的权重为 1。该方法拥有模型灵活、拟合效果好的特点，可以很好的拟合传播数 。其形式如下：

其中 。

为了后续估计参数方便，将 改写为概率密度函数

记 ，则概率密度函数为

其中， 满足约束 , ；；, , 均为估计量。

在混合三参数 logistic 模型中，参数 的估计方法有很多种，其中最大似然估计法由于具有吸引人的渐近统计特性而成为最常用的估计方法之一。期望最大化（Expectation-Maximum，EM）算法是将正态分布模型的混合模型拟合到观测数据的标准方法，其收敛于混合参数的极大似然估计。

Dempster 和 Rubin 在 1977 年提出了 EM 算法²，它是一种通过给定的缺失参数以计算概率密度函数（E 步）和利用计算结果对参数重新进行估计（M 步）交替进行的迭代算法。通过 EM 算法估计混合三参数 logistic 分布模型参数（即 , , ）的步骤如下：

1. 定义 值，在本研究中，；另外设置 , , 的初始值分别为 , , 。

2. E 步 —— 计算似然函数

   记样本序列 共有 个样本，则有

   取对数有

3. M 步 —— 计算后验概率和 , ,

方便起见，先计算 ，即 ：

由于 满足约束 ，故使用 Lagrange 乘数法进行计算。构造 Lagrange 函数

令

记 ，这里 即为 Bayesian 后验概率，则

显然 。另外考虑到 ，则

解得 ，带回原式得

接下来计算 ，即 ：

令

综上

4. 重复 E 步和 M 步，直到参数收敛到所需精度。

我们使用了 python 3.9.13³ 和 scipy 1.7.3⁴ 以应用上述算法。

EM 算法的评估

均方根误差（RMSE）和相关系数（）用于评估EM算法的混合logistic模型的误差。当RMSE越小及相关系数越接近1时结果越好：

其中 为样本值向量，其中 表示第 个样本的样本值； 为估计值向量，其中 表示第 i 个样本的估计值； 表示随机向量 的协方差， 表示随机向量 的方差。

-means 时间序列聚类

-means 是一种基于数据和聚类质心之间的距离度量对数据进行的算法，这种方法对于时间序列数据同样可用。⁵⁶相比于其他算法而言，-means 算法逻辑简单、效率较高、收敛速度快。

考虑到在不同地区，SARS-CoV-2和alpha变种的相对传播能力不同，我们使用了基于动态时间修正（Dynamic Time Warping，DTW）方法⁷的k-means聚类算法对英格兰的所有地区、苏格兰、威尔士和北爱尔兰从2020-03-31到2021-02-28这段时间的SARS-CoV-2和alpha变种的相对传播能力进行了考察，方法如下：

算法如下所示：

1. 在所有样本中随机选取 个地区分别作为 个聚类的中心，在本研究中，。

2. 分配步：将所有样本按其与其最接近（通过 DTW 距离度量）的聚类中心按照数量均等分配给 个群组，即根据这种方法生成的 Voronoi 图对所有样本进行划分，即

   其中 表示第 个聚类， 表示第 和 个聚类的中心， 表示时间序列向量 与 的 DTW 距离，其定义如下：

   其中 表示 和 的 Euclidean 距离度量， 为一条满足以下条件的路径：

   • ,
   • ,
   • .

3. 更新步：根据已分配的样本，重新计算每个聚类的中心，即

4. 重复分配步和更新步，直到所有聚类元素不再变动为止。

上述计算和模型建立过程通过 python 3.9.13³ 和 tslearn 0.5.2⁸ 包完成。

地图

所有的地图均来自来自于英国国家统计署⁹，并通过 Open Government Licence v3.0 协议辅助 GIS 软件 QGIS 3.22¹⁰ 绘制，其中地图投影方法为横轴墨卡托投影¹¹，坐标系统为 EPSG:27700 - OSGB 1936 / British National Grid。

参考文献

Footnotes

1. Birnbaum, A. L. (1968). Some latent trait models and their use in inferring an examinee's ability. Statistical theories of mental test scores. ↩

2. Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1), 1-22. ↩

3. Van Rossum, G., & Drake Jr, F. L. (1995). Python tutorial (Vol. 620). Amsterdam, The Netherlands: Centrum voor Wiskunde en Informatica. ↩ ↩²

4. Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., ... & Van Mulbregt, P. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature methods, 17(3), 261-272. ↩

5. Hartigan, J. A., & Wong, M. A. (1979). Algorithm AS 136: A k-means clustering algorithm. Journal of the royal statistical society. series c (applied statistics), 28(1), 100-108. ↩

6. Huang, X., Ye, Y., Xiong, L., Lau, R. Y., Jiang, N., & Wang, S. (2016). Time series k-means: A new k-means type smooth subspace clustering for time series data. Information Sciences, 367, 1-13. ↩

7. Sakoe, H., & Chiba, S. (1978). Dynamic programming algorithm optimization for spoken word recognition. IEEE transactions on acoustics, speech, and signal processing, 26(1), 43-49. ↩

8. Tavenard, R., Faouzi, J., Vandewiele, G., Divo, F., Androz, G., Holtz, C., ... & Woods, E. (2020). Tslearn, a machine learning toolkit for time series data. J. Mach. Learn. Res., 21(118), 1-6. ↩

9. Office for National Statistics (2020). Local Authority Districts (December 2020) UK. https://geoportal.statistics.gov.uk/maps/bc2820b03de244698c0b0771ae4f345f. ↩

10. Moyroud, N., & Portet, F. (2018). Introduction to QGIS. QGIS and generic tools, 1, 1-17. ↩

11. Jackson, J. E. (1978). Transverse Mercator projection. Survey Review, 24(188), 278-285. ↩