本研究获得北京航空航天大学第三十一届“冯如杯”学生科技竞赛数理类二等奖。
Birnbaum1
L ( x ; a , b , c ) = a 1 + exp ( b + c x ) + c o n s t L(x; a, b, c)=\frac{a}{1+\exp(b+cx)}+\mathrm{const} L ( x ; a , b , c ) = 1 + exp ( b + c x ) a + const 其导数为：
L ′ ( x ; a , b , c ) = − a c exp ( b + c x ) ( 1 + exp ( b + c x ) ) 2 L'(x; a, b, c)=-\frac{ac\exp(b+cx)}{(1+\exp(b+cx))^2} L ′ ( x ; a , b , c ) = − ( 1 + exp ( b + c x ) ) 2 a c exp ( b + c x ) 混合 三参数 logistic 模型是两个 logistic 模型的线性叠加，且每个模型的权重为 1 。该方法拥有模型灵活、拟合效果好的特点，可以很好的拟合传播数 R 0 R_0 R 0
L ( x ; a , b , c ) = ∑ k = 1 K f ( x ; a k , b k , c k ) L(x; a, b, c)=\sum_{k=1}^Kf(x; a_k, b_k, c_k) L ( x ; a , b , c ) = k = 1 ∑ K f ( x ; a k , b k , c k ) 其中 f ( x ; a k , b k , c k ) = a k exp ( b k + c k x ) ( 1 + exp ( b k + c k x ) ) 2 f(x; a_k, b_k, c_k)=\frac{a_k\exp(b_k+c_kx)}{(1+\exp(b_k+c_kx))^2} f ( x ; a k , b k , c k ) = ( 1 + e x p ( b k + c k x ) ) 2 a k e x p ( b k + c k x )
为了后续估计参数方便，将 L ( x ; a , b , c ) L(x; a, b, c) L ( x ; a , b , c )
p ( x ; λ , a , b ) = ∑ k = 1 K λ k f ( x ; a k , b k ) ∫ − ∞ + ∞ ∑ j = 1 K λ i f ( x ; a j , b j ) = ∑ k = 1 K λ k π ∑ j = 1 K λ j f ( x ; a k , b k ) p(x; \boldsymbol{\lambda}, \boldsymbol{a}, \boldsymbol{b})
=\sum_{k=1}^{K}\frac{\lambda_k f(x;a_k,b_k)}{\int_{-\infty}^{+\infty}\sum_{j=1}^K\lambda_i f(x;a_j,b_j)}
=\sum_{k=1}^{K}\frac{\lambda_k}{\pi\sum_{j=1}^K \lambda_j}f(x;a_k,b_k) p ( x ; λ , a , b ) = k = 1 ∑ K ∫ − ∞ + ∞ ∑ j = 1 K λ i f ( x ; a j , b j ) λ k f ( x ; a k , b k ) = k = 1 ∑ K π ∑ j = 1 K λ j λ k f ( x ; a k , b k ) 记 μ k = λ k π ∑ j = 1 K λ j \mu_k=\frac{\lambda_k}{\pi\sum_{j=1}^K \lambda_j} μ k = π ∑ j = 1 K λ j λ k
p ( x ; μ , a , b ) = ∑ k = 1 K μ k f ( x ; a k , b k ) p(x; \boldsymbol{\mu}, \boldsymbol{a}, \boldsymbol{b})
=\sum_{k=1}^{K}\mu_kf(x;a_k,b_k) p ( x ; μ , a , b ) = k = 1 ∑ K μ k f ( x ; a k , b k ) 其中，μ k \mu_k μ k ∑ k = 1 K μ k = 1 π \sum_{k=1}^{K}\mu_k=\frac{1}{\pi} ∑ k = 1 K μ k = π 1 f ( x ; a k , b k ) = a k exp ( b k + a k x ) ( 1 + exp ( b k + a k x ) ) 2 f(x; a_k, b_k)=\frac{a_k\exp(b_k+a_k x)}{(1+\exp(b_k+a_k x))^2} f ( x ; a k , b k ) = ( 1 + e x p ( b k + a k x ) ) 2 a k e x p ( b k + a k x ) ∫ − ∞ + ∞ p ( x ; μ , a , b ) = 1 \int_{-\infty}^{+\infty}p(x; \boldsymbol{\mu}, \boldsymbol{a}, \boldsymbol{b})=1 ∫ − ∞ + ∞ p ( x ; μ , a , b ) = 1 μ = ( μ 1 ⋯ μ K ) T \boldsymbol{\mu}=(\mu_1\cdots \mu_K)^{\mathrm{T}} μ = ( μ 1 ⋯ μ K ) T a = ( a 1 ⋯ a K ) T \boldsymbol{a}=(a_1\cdots a_K)^{\mathrm{T}} a = ( a 1 ⋯ a K ) T b = ( b 1 ⋯ b K ) T \boldsymbol{b}=(b_1\cdots b_K)^{\mathrm{T}} b = ( b 1 ⋯ b K ) T
在混合三参数 logistic 模型中，参数 μ , a , b \boldsymbol{\mu}, \boldsymbol{a}, \boldsymbol{b} μ , a , b E xpectation-M aximum，EM ）算法是将正态分布模型的混合模型拟合到观测数据的标准方法，其收敛于混合参数的极大似然估计。
Dempster 和 Rubin 在 1977 年提出了 EM 算法2 μ \boldsymbol{\mu} μ a \boldsymbol{a} a b \boldsymbol{b} b
定义 K K K K = 2 K=2 K = 2 μ \boldsymbol{\mu} μ a \boldsymbol{a} a b \boldsymbol{b} b μ { 0 } \boldsymbol{\mu}^{\{0\}} μ { 0 } a { 0 } \boldsymbol{a}^{\{0\}} a { 0 } b { 0 } \boldsymbol{b}^{\{0\}} b { 0 }
E 步 —— 计算似然函数
记样本序列 X = ( x 1 , ⋯ , x n ) T \boldsymbol{X}=(x_1, \cdots, x_n)^{\mathrm{T}} X = ( x 1 , ⋯ , x n ) T n n n
P ( X ; μ { n } , a { n } , b { n } ) = ∏ i = 1 N ( ∑ k = 1 K μ k { n } f ( x i ; a k { n } , b k { n } ) ) P(\boldsymbol{X};\boldsymbol{\mu}^{\{n\}}, \boldsymbol{a}^{\{n\}}, \boldsymbol{b}^{\{n\}})
=\prod_{i=1}^N(\sum_{k=1}^{K}\mu_k^{\{n\}}f(x_i;a_k^{\{n\}},b_k^{\{n\}})) P ( X ; μ { n } , a { n } , b { n } ) = i = 1 ∏ N ( k = 1 ∑ K μ k { n } f ( x i ; a k { n } , b k { n } )) 取对数有
ln P ( X ; μ { n } , a { n } , b { n } ) = ∑ i = 1 N ln ∑ k = 1 K μ k { n } f ( x i ; a k { n } , b k { n } ) \ln P(\boldsymbol{X};\boldsymbol{\mu}^{\{n\}}, \boldsymbol{a}^{\{n\}}, \boldsymbol{b}^{\{n\}})
=\sum_{i=1}^N \ln \sum_{k=1}^{K}\mu_k^{\{n\}}f(x_i;a_k^{\{n\}},b_k^{\{n\}}) ln P ( X ; μ { n } , a { n } , b { n } ) = i = 1 ∑ N ln k = 1 ∑ K μ k { n } f ( x i ; a k { n } , b k { n } )
M 步 —— 计算后验概率和 μ { n + 1 } \boldsymbol{\mu}^{\{n+1\}} μ { n + 1 } a { n + 1 } \boldsymbol{a}^{\{n+1\}} a { n + 1 } b { n + 1 } \boldsymbol{b}^{\{n+1\}} b { n + 1 }
方便起见，先计算 μ { n + 1 } \boldsymbol{\mu}^{\{n+1\}} μ { n + 1 } μ k { n + 1 } \mu_k^{\{n+1\}} μ k { n + 1 }
由于 μ k \mu_k μ k ∑ k = 1 K μ k = 1 π \sum_{k=1}^{K}\mu_k=\frac{1}{\pi} ∑ k = 1 K μ k = π 1
L ( X ; μ , a , b ; σ ) = ∑ i = 1 N ln ∑ k = 1 K μ k f ( x i ; a k , b k ) − σ ( ∑ k = 1 K μ k − 1 π ) \mathcal{L}(\boldsymbol{X};\boldsymbol{\mu}, \boldsymbol{a}, \boldsymbol{b};\sigma)
=\sum_{i=1}^N \ln \sum_{k=1}^{K}\mu_kf(x_i;a_k,b_k)-\sigma(\sum_{k=1}^{K}\mu_k-\frac{1}{\pi}) L ( X ; μ , a , b ; σ ) = i = 1 ∑ N ln k = 1 ∑ K μ k f ( x i ; a k , b k ) − σ ( k = 1 ∑ K μ k − π 1 ) 令
∂ L ∂ μ k = ∑ i = 1 N f ( x i ; a k , b k ) ∑ j = 1 K μ j f ( x i ; a j , b j ) − σ = 0 \frac{\partial \mathcal{L}}{\partial \mu_k}
=\sum_{i=1}^N \frac{f(x_i;a_k,b_k)}{\sum_{j=1}^{K}\mu_j f(x_i;a_j,b_j)}-\sigma
=0 ∂ μ k ∂ L = i = 1 ∑ N ∑ j = 1 K μ j f ( x i ; a j , b j ) f ( x i ; a k , b k ) − σ = 0 记 γ i k = μ k f ( x i ; a k , b k ) ∑ j = 1 K μ j f ( x i ; a j , b j ) \gamma_{ik}=\frac{\mu_kf(x_i;a_k,b_k)}{\sum_{j=1}^{K}\mu_j f(x_i;a_j,b_j)} γ ik = ∑ j = 1 K μ j f ( x i ; a j , b j ) μ k f ( x i ; a k , b k ) γ i k \gamma_{ik} γ ik
σ = ∑ i = 1 N γ i k μ k \sigma=\frac{\sum_{i=1}^N\gamma_{ik}}{\mu_k} σ = μ k ∑ i = 1 N γ ik 显然 μ k ∝ ∑ i = 1 n γ i k \mu_k \propto \sum_{i=1}^n\gamma_{ik} μ k ∝ ∑ i = 1 n γ ik ∑ k = 1 K μ k = 1 π \sum_{k=1}^{K}\mu_k=\frac{1}{\pi} ∑ k = 1 K μ k = π 1
∑ k = 1 K μ k = ∑ k = 1 K ∑ i = 1 N γ i k σ = 1 σ ∑ k = 1 K ∑ i = 1 N γ i k = 1 σ ∑ i = 1 N ∑ k = 1 K γ i k = 1 σ ∑ i = 1 N ∑ k = 1 K μ k f ( x i ; a k , b k ) ∑ j = 1 K μ j f ( x i ; a j , b j ) = N σ = 1 π \sum_{k=1}^{K}\mu_k
=\sum_{k=1}^{K}\frac{\sum_{i=1}^N\gamma_{ik}}{\sigma}
=\frac{1}{\sigma}\sum_{k=1}^{K}\sum_{i=1}^N\gamma_{ik}
=\frac{1}{\sigma}\sum_{i=1}^N\sum_{k=1}^{K}\gamma_{ik}
=\frac{1}{\sigma}\sum_{i=1}^N\frac{\sum_{k=1}^{K}\mu_kf(x_i;a_k,b_k)}{\sum_{j=1}^{K}\mu_j f(x_i;a_j,b_j)}
=\frac{N}{\sigma}
=\frac{1}{\pi} k = 1 ∑ K μ k = k = 1 ∑ K σ ∑ i = 1 N γ ik = σ 1 k = 1 ∑ K i = 1 ∑ N γ ik = σ 1 i = 1 ∑ N k = 1 ∑ K γ ik = σ 1 i = 1 ∑ N ∑ j = 1 K μ j f ( x i ; a j , b j ) ∑ k = 1 K μ k f ( x i ; a k , b k ) = σ N = π 1 解得 σ = N π \sigma=N\pi σ = N π
μ k = 1 N π ∑ i = 1 N γ i k \mu_k=\frac{1}{N\pi}\sum_{i=1}^N\gamma_{ik} μ k = N π 1 i = 1 ∑ N γ ik 接下来计算 a { n + 1 } \boldsymbol{a}^{\{n+1\}} a { n + 1 } a k { n + 1 } a_k^{\{n+1\}} a k { n + 1 }
令
∂ ln P ∂ a k = \frac{\partial \ln P}{\partial a_k}
= ∂ a k ∂ ln P = 综上
\displaylines { γ i k { n } = μ k { n } f ( x i ; a k { n } , b k { n } ) ∑ j = 1 K μ j { n } f ( x i ; a j { n } , b j { n } ) μ k { n + 1 } = 1 N π ∑ i = 1 N γ i k { n } a k { n + 1 } b k { n + 1 } \displaylines{
\left\{
\begin{array}{lr}
\displaystyle \gamma_{ik}^{\{n\}}=\frac{\mu_k^{\{n\}}f(x_i;a_k^{\{n\}},b_k^{\{n\}})}{\sum_{j=1}^{K}\mu_j^{\{n\}} f(x_i;a_j^{\{n\}},b_j^{\{n\}})} \\
\displaystyle \mu_k^{\{n+1\}}=\frac{1}{N\pi}\sum_{i=1}^N\gamma_{ik}^{\{n\}} \\
\displaystyle a_k^{\{n+1\}} \\
\displaystyle b_k^{\{n+1\}}
\end{array}
\right.
} \displaylines ⎩ ⎨ ⎧ γ ik { n } = ∑ j = 1 K μ j { n } f ( x i ; a j { n } , b j { n } ) μ k { n } f ( x i ; a k { n } , b k { n } ) μ k { n + 1 } = N π 1 i = 1 ∑ N γ ik { n } a k { n + 1 } b k { n + 1 }
重复 E 步 和 M 步 ，直到参数收敛到所需精度。
我们使用了 python 3.9.133 4
均方根误差（RMSE）和相关系数（ρ \rho ρ
\displaylines e i = y i ^ − y i R M S E = 1 n ∑ i = 1 n e i 2 ρ = c o v ( y ^ , y ) v a r ( y ^ ) ⋅ v a r ( y ) \displaylines{
e_i=\hat{y_i}-y_i \\
\mathrm{RMSE}=\sqrt{\frac{1}{n}\sum_{i=1}^ne_i^2} \\
\rho=\frac{\mathrm{cov}(\hat{y}, y)}{\sqrt{\mathrm{var}(\hat{y})\cdot\mathrm{var}(y)}}
} \displaylines e i = y i ^ − y i RMSE = n 1 i = 1 ∑ n e i 2 ρ = var ( y ^ ) ⋅ var ( y ) cov ( y ^ , y ) 其中 y = ( y 1 , y 2 , ⋯ , y i , ⋯ , y n ) T y=(y_1,y_2,\cdots,y_i,\cdots,y_n)^T y = ( y 1 , y 2 , ⋯ , y i , ⋯ , y n ) T y i y_i y i i i i y ^ = ( y ^ 1 , y ^ 2 , ⋯ , y ^ i , ⋯ , y ^ n ) T \hat{y}=(\hat{y}_1,\hat{y}_2,\cdots,\hat{y}_i,\cdots,\hat{y}_n)^T y ^ = ( y ^ 1 , y ^ 2 , ⋯ , y ^ i , ⋯ , y ^ n ) T y ^ i \hat{y}_i y ^ i c o v ( a , b ) \mathrm{cov}(a,b) cov ( a , b ) a , b a,b a , b v a r ( a ) var(a) v a r ( a ) a a a
k k k 5 6 k k k
考虑到在不同地区，SARS-CoV-2和alpha变种的相对传播能力不同，我们使用了基于动态时间修正（Dynamic Time Warping，DTW）方法7
算法如下所示：
在所有样本中随机选取 K K K K K K K = 2 K=2 K = 2
分配步 ：将所有样本按其与其最接近（通过 DTW 距离度量）的聚类中心按照数量均等分配给 K K K
S i = { x p : D T W ( x p , m i ) ≤ D T W ( x p , m j ) ∀ 1 ≤ j ≤ K } S_i=\{\boldsymbol{x}_p: \mathrm{DTW}(\boldsymbol{x}_p, \boldsymbol{m}_i)\leq\mathrm{DTW}(\boldsymbol{x}_p, \boldsymbol{m}_j)\qquad \forall\,1\leq j\leq K\} S i = { x p : DTW ( x p , m i ) ≤ DTW ( x p , m j ) ∀ 1 ≤ j ≤ K } 其中 S i S_i S i i i i m i , m j \boldsymbol{m}_i, \boldsymbol{m}_j m i , m j i i i j j j D T W ( a , b ) \mathrm{DTW}(\boldsymbol{a}, \boldsymbol{b}) DTW ( a , b ) a = ( a 1 , ⋯ , a n ) \boldsymbol{a}=(a_1,\cdots,a_n) a = ( a 1 , ⋯ , a n ) b = ( b 1 , ⋯ , b m ) \boldsymbol{b}=(b_1,\cdots,b_m) b = ( b 1 , ⋯ , b m )
D T W ( a , b ) = min π ∑ ( i , j ) ∈ π ∥ a i − b j ∥ 2 \mathrm{DTW}(\boldsymbol{a}, \boldsymbol{b})=\min_{\boldsymbol{\pi}}\sqrt{\sum_{(i,j)\in\boldsymbol{\pi}}\Vert a_i-b_j \Vert^2} DTW ( a , b ) = π min ( i , j ) ∈ π ∑ ∥ a i − b j ∥ 2 其中 ∥ a i − b j ∥ \Vert a_i-b_j \Vert ∥ a i − b j ∥ a i a_i a i b j b_j b j π = { π 0 , ⋯ , π S } \boldsymbol{\pi}=\{\pi_0,\cdots,\pi_S\} π = { π 0 , ⋯ , π S }
∀ π s = ( i s , j s ) , 0 ≤ i s ≤ n , 0 ≤ j s ≤ m \forall\,\pi_s=(i_s, j_s), 0\leq i_s\leq n, 0\leq j_s\leq m ∀ π s = ( i s , j s ) , 0 ≤ i s ≤ n , 0 ≤ j s ≤ m π 0 = ( 0 , 0 ) , π S = ( n − 1 , m − 1 ) \pi_0=(0,0), \pi_S=(n-1,m-1) π 0 = ( 0 , 0 ) , π S = ( n − 1 , m − 1 ) i k − 1 ≤ i k ≤ i k + 1 , j k − 1 ≤ j k ≤ j k + 1 i_{k-1}\leq i_k\leq i_{k+1}, j_{k-1}\leq j_k\leq j_{k+1} i k − 1 ≤ i k ≤ i k + 1 , j k − 1 ≤ j k ≤ j k + 1
更新步 ：根据已分配的样本，重新计算每个聚类的中心，即
m i = 1 ∣ S i ∣ ∑ x j ∈ S i x j \boldsymbol{m}_i=\frac{1}{|S_i|}\sum_{\boldsymbol{x}_j\in S_i}\boldsymbol{x}_j m i = ∣ S i ∣ 1 x j ∈ S i ∑ x j
重复分配步和更新步，直到所有聚类元素不再变动为止。
上述计算和模型建立过程通过 python 3.9.133 8
所有的地图均来自来自于英国国家统计署9 10 11