零一膨胀几何分布的统计分析及应用

刘梦瑶; 肖翔

doi:10.12299/jsues.23-0213

摘要: 研究了0−1膨胀几何分布模型，构造隐变量的条件分布，并设计抽样算法。在数据扩充的基础上，运用极大似然估计、期望极大（expectation maximization，EM）算法及贝叶斯方法对模型参数进行估计。设定不同的样本量和参数真值，通过数值模拟对上述方法进行性能评估。最后，对1994年美国底特律交通事故死亡数据集进行分析，研究表明，0−1膨胀几何分布模型能够较好地对该数据集进行拟合。

关键词:

Abstract: Zero-and-one-inflated geometric distribution model was investigated, the conditional distribution of latent variables were constructed, and a sampling algorithm was designed. On the basis of data expansion, the maximum likelihood estimation, expectation maximization (EM) algorithm and Bayesian method were employed to estimate the model parameters. Different sample sizes and parameter true values were setted, and the performance of these methods were evaluated through numerical simulations. Finally, the Detroit traffic accident deaths dataset in 1994 of United States were analyzed, the results indicate that zero-and-one-inflated geometric distribution model can fit the dataset better.

Key words:

zero-and-one-inflated geometric distribution /
data augmentation /
maximum likelihood estimation /
expectation maximization (EM) algorithm /
Bayesian estimation

计数数据是统计学研究的热点，具有过多零值的过度分散计数数据经常出现在自然科学、医学以及保险等各个领域。尽管泊松模型被广泛用于分析此类离散数据，但其均值和方差相等的强假设意味着当数据的方差大于均值时，该模型不再适用于对这类数据进行建模，这种现象被称为“过度分散”。

过度分散的一个常见来源是观察到的计数包含过多的零值，因此，零膨胀模型成为众多学者研究的对象。张良超等^[1]通过零膨胀泊松模型讨论了3种损失函数下风险参数的贝叶斯估计。肖翔等^[2]提出零膨胀几何分布模型并进行参数估计。但样本数据有时会出现大量的零值和一值。例如，电子商务已深入百姓生活，人们在大型商场购物时往往抱着货比三家的心态，很多顾客选择不购买衣服或者只购买一件衣服。对于这类数据集，零膨胀模型不适合进行建模，需要进一步将其扩展到零一膨胀模型，常见的模型有：零一膨胀泊松（ZOIP）模型和零一膨胀负二项（ZOINB）模型。

Zhang等^[3]构建了ZOIP模型的5个等价表达式，并推导出重要的分布性质。Tang等^[4]建立了ZOIP模型的新型结构，通过期望极大（expectation maximization，EM）算法和贝叶斯方法进行参数估计，数值模拟的效果良好。Jorinsian等^[5]提出零一膨胀负二项−贝塔（ZOINB-BE）分布并研究其重要性质，通过该模型分析了3个真实数据集。Tajuddin等^[6]介绍了零一膨胀泊松−林德利（ZOIP-Lindley）分布，利用矩估计(moment estimation, MM)和最大似然估计(maximum likelihood estimation, MLE)进行参数估计。研究表明，ZOIP-Lindley分布模型能够较好地描述新西兰白兔因分娩死亡的数量。

目前，对0−1膨胀几何分布及回归模型的研究较少。Xiao等^[7]构建了零一膨胀几何分布（ZOIGE）回归模型，在贝叶斯推断中引入了Pólya-gamma潜变量，将高维参数抽样问题转化为低维参数抽样和潜变量抽样，实例分析表明其比ZOIP回归模型更适合分析博士发表论文数量的数据集。为了更好地拟合零一值过多的数据，通过重参数化，肖翔等^[8]对0−1膨胀几何分布模型进行客观贝叶斯分析，得到Jeffreys先验以及reference先验，达到很好的拟合效果，然后设计抽样算法对后验分布进行抽样，获得参数的估计。本研究以0−1膨胀几何分布模型为研究对象，分别使用极大似然估计、EM算法下的极大似然估计和贝叶斯估计3种方法进行参数估计，通过数值仿真和实例分析，评价这3种方法的估计效果，以提高鉴别真实零、一观测值的效率，丰富零一膨胀模型的研究方法。

1. 0−1膨胀几何分布

在伯努利试验中，假设每次试验成功的概率为 $\theta$ ， $Y$ 为直到第一次成功出现为止，所进行的试验次数，其概率分布为 $P\{ Y = k\} = {(1 - \theta )^{k - 1}}\theta$ ， $k = 1, 2, \cdots$ ，称 $Y$ 服从几何分布。

换个角度进行定义， $Y$ 为直到第一次失败出现为止，试验成功的次数，其概率分布为 $P\{ Y = k\} = {\theta ^k}(1 - \theta )$ , $k = 0,1, \cdots$ ，称 $Y$ 服从几何分布。这两种定义是等价的，本研究采用第2种几何分布定义。

当样本数据出现大量的零值和一值，且数据的尾部退化较快时，为了更好地描述这种现象，本研究构建零一膨胀几何分布模型，一个非负ZOIGE分布变量 $Z = B(1 - X) + (1 - B)Y$ 。其中， $(B,X,Y)$ 相互独立， $B$ 服从于试验成功概率为 $p$ 的伯努利分布， $X$ 服从于试验成功概率为 $q$ 的伯努利分布， $Y$ 服从试验成功概率为 $\theta$ 的几何分布，即 $P\{ Y = k\} = {\theta ^k}(1 - \theta )$ （ $k = 0,1 \cdots$ ）。 $(B,X,Y)$ 无法直接观测，称为隐变量， $Z$ 与 $(B,X,Y)$ 的关系为

$\left\{ \begin{gathered} \{ Z = 0\} \Leftrightarrow \{ Y = 0,B = 0\} \cup \{ X = 1,B = 1\} \\ \{ Z = 1\} \Leftrightarrow \{ Y = 1,B = 0\} \cup \{ X = 0,B = 1\} \\ \{ Z = k\} \Leftrightarrow \{ Y = k,B = 0\} ,\;\; k = 2,3 \cdots \\ \end{gathered} \right.$

(1)

根据式(1)，随机变量 $Z$ 的概率分布为

$P\{ Z = k\} = \left\{ \begin{gathered} pq + (1 - p)(1 - \theta ), \;\; k = 0 \\ p(1 - q) + (1 - p)\theta (1 - \theta ), \;\; k = 1 \\ (1 - p){\theta ^k}(1 - \theta ), \;\; k = 2,3, \cdots \\ \end{gathered} \right.$

(2)

式中： $0 \leqslant p \leqslant 1$ ， $0 \leqslant q \leqslant 1$ ， $0 \leqslant \theta \leqslant 1$ 。根据式(2)， $p$ 和 $\left( {1 - p} \right)$ 分别为伯努利分布和几何分布的混合比例。当 $p = 0$ 时，模型为几何分布；当 $q = 1$ 时，模型为零膨胀几何分布(ZIGE)模型。

2. 极大似然估计

设 $Z = ({Z_1},{Z_2}, \cdots ,{Z_n})$ 为ZOIGE分布的样本观测值，由式(2)可得， $(p,q,\theta )$ 的似然函数为

$\begin{split} & L(p,q,\theta |Z) = {[pq + (1 - p)(1 - \theta )]^{{M_0}}} \cdot [p(1 - q) + \\ & \quad (1 - p)\theta (1 - \theta )]^{{M_1}} \cdot {(1 - p)^{n - {M_0} - {M_1}}}{\theta ^M}\cdot \\ & \quad (1 -\theta )^{n - {M_0} - {M_1}} \end{split}$

(3)

式中： ${M_0}$ 为集合 $\{ i:{Z_i} = 0\}$ 的元素个数； ${M_1}$ 为集合 $\{ i:{Z_i} = 1\}$ 的元素个数； $M =\displaystyle \sum\nolimits_{{Z_i} \geqslant 2} {{Z_i}}$ 。式(3)中，令

$\left\{ \begin{gathered} \varphi = pq + (1 - p)(1 - \theta ) \\ \phi = p(1 - q) + (1 - p)\theta (1 - \theta ) \\ \end{gathered} \right.$

(4)

由式(4)可得， $1 - p = \dfrac{{1 - \varphi - \phi }}{{{\theta ^2}}}$ 。经过重参数化后， $(\varphi ,\phi , \theta )$ 的似然函数为

$\begin{split} & L(\varphi ,\phi ,\theta |Z) = \\ & \quad {\varphi ^{{M_0}}}{\phi ^{{M_1}}}{(1 - \varphi - \phi )^{n - {M_0} - {M_1}}}\dfrac{{{\theta ^M}{{(1 - \theta )}^{n - {M_0} - {M_1}}}}}{{{\theta ^{2(n - {M_0} - {M_1})}}}} \end{split}$

分别令 $\dfrac{{\partial \ln L(\varphi ,\phi ,\theta |Z)}}{{\partial \varphi }} = 0\text{，}\dfrac{{\partial \ln L(\varphi ,\phi ,\theta |Z)}}{{\partial \phi }} = 0\text{，} \dfrac{{\partial \ln L(\varphi ,\phi ,\theta |Z)}}{{\partial \theta }} = 0$ ，得到 $(\varphi ,\phi ,\theta )$ 的极大似然估计为

$\hat \varphi = \dfrac{{{M_0}}}{n}\text{，}\hat \phi = \dfrac{{{M_1}}}{n}\text{，}\hat \theta = \dfrac{{M - 2(n - {M_0} - {M_1})}}{{M - (n - {M_0} - {M_1})}}。$

根据极大似然估计的不变性，得到 $p,q$ 的极大似然估计为

$\hat p = \dfrac{{\hat \varphi + \hat \phi + {{\hat \theta }^2} - 1}}{{{{\hat \theta }^2}}}\text{，}\hat q = \dfrac{{\hat \varphi - (1 - \hat p)(1 - \hat \theta )}}{{\hat p}}。$

3. 潜变量的分布及抽样

3.1 潜变量的条件概率分布

由式(1)推导潜变量 $(B,X,Y)$ 的条件概率分布为 $\pi (B,X,Y|Z,p,q,\theta )$ 。

当 $B = 0，X = 0，Y = 0$ 时，由 $B = 0，Y = 0$ 可以推出 $Z = 0$ ，则

$\begin{split} & P\{ B = 0, X = 0, Y = 0|Z = 0, p, q, \theta \} = \\ &\quad\dfrac{{P\{ B = 0, X = 0, Y = 0, Z = 0\} }}{{P\{ Z = 0\} }} =\\ &\quad \dfrac{{P\{ B = 0, X = 0, Y = 0\} }}{{P\{ Z = 0\} }} = \\ &\quad\dfrac{{P\{ B = 0\} \cdot P\{ X = 0\} \cdot P\{ Y = 0\} }}{{P\{ Z = 0\} }} = \\ &\quad\dfrac{{(1 - p)(1 - q)(1 - \theta )}}{{pq + (1 - p)(1 - \theta )}} \end{split}$

当 $B = 0，X = 1，Y = 0$ 时，由 $B = 0，Y = 0$ 可以推出 $Z = 0$ ，则

$\begin{split} & P\{ B = 0, X = 1, Y = 0|Z = 0, p, q, \theta \} =\\ & \quad \dfrac{{P\{ B = 0, X = 1, Y = 0, Z = 0\} }}{{P\{ Z = 0\} }} =\\ & \quad \dfrac{{P\{ B = 0, X = 1, Y = 0\} }}{{P\{ Z = 0\} }} = \\ & \quad \dfrac{{P\{ B = 0\} \cdot P\{ X = 1\} \cdot P\{ Y = 0\} }}{{P\{ Z = 0\} }} = \\ & \quad \dfrac{{(1 - p)q(1 - \theta )}}{{pq + (1 - p)(1 - \theta )}} \end{split}$

当 $B = 1，X = 1，Y = k$ $(k = 0,1 \cdots )$ 时，由 $B = 1， X = 1$ 可以推出 $Z = 0$ ，则

$\begin{split} & P\{ B = 1,X = 1,Y = k|Z = 0,p,q,\theta \} = \\ & \quad \dfrac{{P\{ B = 1,X = 1,Y = k,Z = 0\} }}{{P\{ Z = 0\} }} =\\ & \quad \dfrac{{P\{ B = 1,X = 1,Y = k\} }}{{P\{ Z = 0\} }} = \\ & \quad \dfrac{{P\{ B = 1\} \cdot P\{ X = 1\} \cdot P\{ Y = k\} }}{{P\{ Z = 0\} }} = \\ & \quad \dfrac{{pq{\theta ^k}(1 - \theta )}}{{pq + (1 - p)(1 - \theta )}} \\ \end{split}$

易证以上3个条件概率之和为 $1$ ，因此有

$\begin{split} & P\{B = i,X = j,Y = k | Z = 0,p,q,\theta \} = \\ & \quad \left\{\begin{array}{ll}\dfrac{(1 - p)(1 - q)(1 - \theta )}{pq + (1 - p)(1 - \theta )},\;\; i = 0,j = 0,k = 0\\ \dfrac{(1 - p)q(1 - \theta )}{pq + (1 - p)(1 - \theta )},\;\; i = 0,j = 1,k = 0\\ \dfrac{pq{\theta }^{k}(1 - \theta )}{pq + (1 - p)(1 - \theta )},\;\; i = 1,j = 1,k = 0,1\cdots \\ 0,\;\; 其他 \end{array}\right. \end{split}$

(5)

同理可得

$\begin{split} & P\{B = i,X = j,Y = k|Z = 1,p,q,\theta \} = \\ & \quad \left\{\begin{array}{l}\dfrac{(1 - p)(1 - q)\theta (1 - \theta )}{p(1 - q) + (1 - p)\theta (1 - \theta )},\;\; i = 0,j = 0,k = 1\\ \dfrac{(1 - p)q\theta (1 - \theta )}{p(1 - q) + (1 - p)\theta (1 - \theta )},\;\; i = 0,j = 1,k = 1\\ \dfrac{p(1 - q){\theta }^{k}(1 - \theta )}{p(1 - q) + (1 - p)\theta (1 - \theta )},\;\; i = 1,j = 0,k = 0,1\cdots \\ 0,\;\; 其他\end{array}\right. \end{split}$

(6)

$\begin{split} & P\{B=i,X=j,Y=k|Z=m,p,q,\theta \}=\\ & \quad \left\{\begin{array}{cc}1-q,& i=0,j=0,k=m\\ q,& i=0,j=1,k=m\\ 0,& 其他\end{array} \right. \end{split}$

(7)

其中 $m = 2,3 \cdots$ 则称式(5)至式(7)为隐变量 $(B,X,Y)$ 的条件概率分布。

3.2 潜变量的抽样方法

由式(1)和条件概率分布 $\pi (B,X,Y|Z,p,q,\theta )$ ，在给定 ${Z_i}$ 条件下，设计 $(B,X,Y)$ 样本的抽样机制。

1）若 ${Z_i} = 0$ ，抛掷一枚正面朝上概率为 $\dfrac{{pq}}{{pq + (1 - p)(1 - \theta )}}$ 的硬币，当出现正面朝上时，令 ${B_i} = 1$ ， ${X_i} = 1$ ， ${Y_i}$ 由几何分布抽样得到；当出现反面朝上时，令 ${B_i} = 0$ ， ${Y_i} = 0$ 。此时，再抛掷另一枚正面朝上概率为 $q$ 的硬币，若该硬币出现正面朝上，则令 ${X_i} = 1$ ；否则，令 ${X_i} = 0$ 。

2）若 ${Z_i} = 1$ ，抛掷一枚正面朝上概率为 $\dfrac{{p(1 - q)}}{{p(1 - q) + (1 - p)\theta (1 - \theta )}}$ 的硬币，当出现正面朝上时，令 ${B_i} = 1$ ， ${X_i} = 0$ ， ${Y_i}$ 由几何分布抽样得到；当出现反面朝上时，令 ${B_i} = 0$ ， ${Y_i} = 1$ 。此时，再抛掷另一枚正面朝上概率为 $q$ 的硬币，若该硬币出现正面朝上，则令 ${X_i} = 1$ ；否则，令 ${X_i} = 0$ 。

3）若 ${Z_i} = m$ $(m = 2,3 \cdots )$ ，抛掷一枚正面朝上概率为 $q$ 的硬币，当出现正面朝上时，令 ${B_i} = 0$ ， ${X_i} = 1$ ， ${Y_i} = m$ ；当出现反面朝上时，令 ${B_i} = 0$ ， ${X_i} = 0$ ， ${Y_i} = m$ 。

构建条件概率分布 $\pi (B,X,Y|Z,p,q,\theta )$ 及其抽样方法是本研究的创新点和重点，在后续的EM算法和贝叶斯推断中起到关键作用。

4. EM算法下的极大似然估计

设 $Z = ({Z_1},{Z_2}, \cdots ,{Z_n})$ 为来自ZOIGE分布的样本观测值， $B = ({B_1},{B_2}, \cdots ,{B_n})$ ， $X = ({X_1},{X_2}, \cdots ,{X_n})$ ， $Y = ({Y_1},{Y_2}, \cdots ,{Y_n})$ 为隐变量组，则数据扩充下 $(p,q,\theta )$ 的完全似然函数为

$\begin{split} & L(p,q,\theta |B,X,Y,Z) = \\ &\quad \prod\limits_{i = 1}^n {{{[p{q^{{X_i}}}{{(1 - q)}^{1 - {X_i}}}]}^{{B_i}}}{{[(1 - p){\theta ^{{Y_i}}}(1 - \theta )]}^{1 - {B_i}}}} =\\ & \quad \prod\limits_{i = 1}^n {p^{{B_i}}}{{(1 - p)}^{1 - {B_i}}}{q^{{B_i}{X_i}}}{{(1 - q)}^{{B_i}(1 - {X_i})}}\cdot\\&\quad {\theta ^{(1 - {B_i}){Y_i}}}{{(1 - \theta )}^{(1 - {B_i})}} \end{split}$

(8)

完全对数似然函数为

$\begin{split} & l(p,q,\theta |B,X,Y,Z) =\displaystyle \sum\limits_{i = 1}^n {[{B_i}\ln p + (1 - {B_i})\ln (1 - p)]} +\\ & \quad \displaystyle \sum\limits_{i = 1}^n {[{B_i}{X_i}\ln q + {B_i}(1 - {X_i})\ln (1 - q)]} +\\ & \quad \displaystyle \sum\limits_{i = 1}^n {(1 - {B_i})[{Y_i}\ln } \theta + \ln (1 - \theta )] \end{split}$

令

$\left\{ {\begin{array}{*{20}{l}} {l(p|B,X,Y,Z) =\displaystyle \sum\limits_{i = 1}^n {[{B_i}\ln p + (1 - {B_i})\ln (1 - p)]} } \\ {l(q|B,X,Y,Z) =\displaystyle \sum\limits_{i = 1}^n {[{B_i}{X_i}\ln q + {B_i}(1 - {X_i})\ln (1 - q)]} } \\ {l(\theta |B,X,Y,Z) =\displaystyle \sum\limits_{i = 1}^n {(1 - {B_i})[{Y_i}\ln } \theta + \ln (1 - \theta )]} \end{array}} \right.$

在EM算法中，在数据扩充策略下， $(p,q,\theta )$ 的极大似然估计由E步与M步交替迭代得到。在E步中，根据当前的参数估计 $(p,q,\theta )$ ，计算隐变量 $({B_i},{X_i},{Y_i})$ 的期望；在M步中，将E步中得到的 $({B_i},{X_i},{Y_i})$ 期望视为常数，对 $l(p|B,X,Y,Z)$ 、 $l(q|B,X,Y,Z)$ 和 $l(\theta |B,X,Y,Z)$ 进行极大化。根据式(5)至式(7)，EM算法具体实施如下。

1）E步：计算 $B_i^{(t + 1)} = E({B_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}})$

当 ${Z_i} = 0$ 时，有

$\begin{split} & E({B_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}) = \\ &\quad P\{ {B_i} = 1|{Z_i} = 0,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} =\\ &\quad \displaystyle \sum\limits_{k = 0}^\infty {P\{ {B_i} = 1,{X_i} = 1,{Y_i} = k|{Z_i} = 0\} } =\\ &\quad \displaystyle \sum\limits_{k = 0}^\infty {\dfrac{{{p^{(t)}}{q^{(t)}}{\theta ^{(t)}}^k(1 - {\theta ^{(t)}})}}{{{p^{(t)}}{q^{(t)}} + (1 - {p^{(t)}})(1 - {\theta ^{(t)}})}}} =\\ &\quad \dfrac{{{p^{(t)}}{q^{(t)}}}}{{{p^{(t)}}{q^{(t)}} + (1 - {p^{(t)}})(1 - {\theta ^{(t)}})}} \end{split}$

当 ${Z_i} = 1$ 时，有

$\begin{split} & E({B_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}) = \\ &\quad P\{ {B_i} = 1|{Z_i} = 1,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} =\\ &\quad \displaystyle \sum\limits_{k = 0}^\infty {P\{ {B_i} = 1,{X_i} = 0,{Y_i} = k|{Z_i} = 1\} } =\\ &\quad \displaystyle \sum\limits_{k = 0}^\infty {\dfrac{{{p^{(t)}}(1 - {q^{(t)}}){\theta ^{(t)}}^k(1 - {\theta ^{(t)}})}}{{{p^{(t)}}(1 - {q^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}}} = \\ &\quad \dfrac{{{p^{(t)}}(1 - {q^{(t)}})}}{{{p^{(t)}}(1 - {q^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}} \end{split}$

当 ${Z_i} = m$ $(m = 2,3 \cdots )$ 时，有 $E({B_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}) = P\{ {B_i} = 1|{Z_i} = m,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} = 0$ 。 ${B_i}$ 的 $t + 1$ 步估计表达式为

$\begin{split} & B_i^{(t + 1)} =\\ & \left\{ {\begin{array}{*{20}{c}} {\dfrac{{{p^{(t)}}{q^{(t)}}}}{{{p^{(t)}}{q^{(t)}} + (1 - {p^{(t)}})(1 - {\theta ^{(t)}})}},}&{{Z_i} = 0} \\ {\dfrac{{{p^{(t)}}(1 - {q^{(t)}})}}{{{p^{(t)}}(1 - {q^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}},}&{{Z_i} = 1} \\ {0,}&{{Z_i} = 2,3 \cdots } \end{array}} \right. \end{split}$

2）E步：计算 $X_i^{(t + 1)} = E({X_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}})$

当 ${Z_i} = 0$ 时，有

$\begin{split} & E({X_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}) =\\ & \quad P\{ {X_i} = 1|{Z_i} = 0,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} =\\ &\quad \displaystyle \sum\limits_{k = 0}^\infty {P\{ {B_i} = 1,{X_i} = 1,{Y_i} = k|{Z_i} = 0\} } + \\ &\quad P\{ {B_i} = 0,{X_i} = 1,{Y_i} = 0|{Z_i} = 0\} =\\ &\quad \displaystyle \sum\limits_{k = 0}^\infty {\dfrac{{{p^{(t)}}{q^{(t)}}{\theta ^{(t)}}^k(1 - {\theta ^{(t)}})}}{{{p^{(t)}}{q^{(t)}} + (1 - {p^{(t)}})(1 - {\theta ^{(t)}})}}} +\\ & \quad \dfrac{{(1 - {p^{(t)}}){q^{(t)}}(1 - {\theta ^{(t)}})}}{{{p^{(t)}}{q^{(t)}} + (1 - {p^{(t)}})(1 - {\theta ^{(t)}})}} = \end{split}$

$\begin{split} \quad \dfrac{{{p^{(t)}}{q^{(t)}} + (1 - {p^{(t)}}){q^{(t)}}(1 - {\theta ^{(t)}})}}{{{p^{(t)}}{q^{(t)}} + (1 - {p^{(t)}})(1 - {\theta ^{(t)}})}} \end{split}$

当 ${Z_i} = 1$ 时，有

$\begin{split} & E({X_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}) = P\{ {X_i} = 1|{Z_i} = \\ & \quad 1,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\}= P\{ {B_i} = 0,{X_i} = 1,{Y_i} = 1|{Z_i} = 1\} = \\ & \quad \dfrac{{(1 - {p^{(t)}}){q^{(t)}}{\theta ^{(t)}}(1 - {\theta ^{(t)}})}}{{{p^{(t)}}(1 - {q^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}}\end{split}$

当 ${Z_i} = m$ $(m = 2,3 \cdots )$ 时，有

$\begin{split} & E({X_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}) = P\{ {X_i} = 1|{Z_i} = \\ & \quad m,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} = P\{ {B_i} = 0,{X_i} = 1,{Y_i} = m|{Z_i} =\\ & \quad m\} = {q^{(t)}} \end{split}$

${X_i}$ 的 $t + 1$ 步估计表达式为

$\qquad\qquad\qquad\qquad\qquad \begin{split} X_i^{(t + 1)} =\left\{ {\begin{array}{*{20}{c}} {\dfrac{{{p^{(t)}}{q^{(t)}} + (1 - {p^{(t)}}){q^{(t)}}(1 - {\theta ^{(t)}})}}{{{p^{(t)}}{q^{(t)}} + (1 - {p^{(t)}})(1 - {\theta ^{(t)}})}},}&{{Z_i} = 0} \\ {\dfrac{{(1 - {p^{(t)}}){q^{(t)}}{\theta ^{(t)}}(1 - {\theta ^{(t)}})}}{{{p^{(t)}}(1 - {q^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}},}&{{Z_i} = 1} \\ {{q^{(t)}},}&{{Z_i} = 2,3 \cdots } \end{array}} \right. \end{split}$

3）E步：计算 $Y_i^{(t + 1)} = E({Y_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}})$

当 ${Z_i} = 0$ 时，有

$\begin{split} & E({Y_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}) =\\ & \quad \displaystyle \sum\limits_{k = 0}^\infty {kP\{ {Y_i} = k|{Z_i} = 0,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} } = \\ & \quad \displaystyle \sum\limits_{k = 0}^\infty {kP\{ {B_i} = 1,{X_i} = 1,{Y_i} = k|{Z_i} = 0\} } =\\ & \quad \displaystyle \sum\limits_{k = 0}^\infty {\dfrac{{k{p^{(t)}}{q^{(t)}}{\theta ^{(t)}}^k(1 - {\theta ^{(t)}})}}{{{p^{(t)}}{q^{(t)}} + (1 - {p^{(t)}})(1 - {\theta ^{(t)}})}}} =\\ & \quad \dfrac{{{p^{(t)}}{q^{(t)}}{\theta ^{(t)}}}}{{{p^{(t)}}{q^{(t)}}(1 - {\theta ^{(t)}}) + (1 - {p^{(t)}}){{(1 - {\theta ^{(t)}})}^2}}} \end{split}$

当 ${Z_i} = 1$ 时，有

$\begin{split} & E({Y_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}) =\\ & \quad \displaystyle \sum\limits_{k = 0}^\infty {kP\{ {Y_i} = k|{Z_i} = 1,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} } = \\ & \quad P\{ {B_i} = 0,{X_i} = 0,{Y_i} = 1|{Z_i} = 1,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} +\\ & \quad P\{ {B_i} = 0,{X_i} = 1,{Y_i} = 1|{Z_i} = 1,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} +\\ & \quad \displaystyle \sum\limits_{k = 0}^\infty {kP\{ {B_i} = 1,{X_i} = 0,{Y_i} = k|{Z_i} = 1,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} } = \end{split}$

$\begin{split} & \quad \dfrac{{(1 - {p^{(t)}})(1 - {q^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}}{{{p^{(t)}}(1 - {q^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}} +\\ & \quad \dfrac{{(1 - {p^{(t)}}){q^{(t)}}{\theta ^{(t)}}(1 - {\theta ^{(t)}})}}{{{p^{(t)}}(1 - {q^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}} +\\ & \quad \displaystyle \sum\limits_{k = 0}^\infty {\dfrac{{k(1 - {p^{(t)}}){q^{(t)}}{\theta ^{(t)k}}(1 - {\theta ^{(t)}})}}{{{p^{(t)}}(1 - {q^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}}} = \\ & \quad \dfrac{{(1 - {p^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}}{{{p^{(t)}}(1 - {q^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}(1 - {\theta ^{(t)}})}} + \\ & \quad \dfrac{{(1 - {p^{(t)}}){q^{(t)}}{\theta ^{(t)}}}}{{{p^{(t)}}(1 - {q^{(t)}})(1 - {\theta ^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}{{(1 - {\theta ^{(t)}})}^2}}} = \\ & \quad \dfrac{{(1 - {p^{(t)}}){\theta ^{(t)}}{{(1 - {\theta ^{(t)}})}^2} + (1 - {p^{(t)}}){q^{(t)}}{\theta ^{(t)}}}}{{{p^{(t)}}(1 - {q^{(t)}})(1 - {\theta ^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}{{(1 - {\theta ^{(t)}})}^2}}} \end{split}$

当 ${Z_i} = m$ $(m = 2,3 \cdots )$ 时，有

$\begin{split} & E({Y_i}|{Z_i},{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}) =\\ & \quad \displaystyle \sum\limits_{k = 0}^\infty {kP\{ {Y_i} = k|{Z_i} = m,{p^{(t)}},{q^{(t)}},{\theta ^{(t)}}\} } =\\ & \quad mP\{ {B_i} = 0,{X_i} = 0,{Y_i} = m|{Z_i} = m\} + \\ & \quad mP\{ {B_i} = 0,{X_i} = 1,{Y_i} = m|{Z_i} = m\} = \\ & \quad m(1 - {q^{(t)}}) + m{q^{(t)}} = m = {Z_i} \end{split}$

${Y_i}$ 的 $t + 1$ 步估计表达式为

$\qquad\qquad\qquad\qquad\begin{split} Y_i^{(t + 1)} = \left\{ {\begin{array}{*{20}{c}} {\dfrac{{{p^{(t)}}{q^{(t)}}{\theta ^{(t)}}}}{{{p^{(t)}}{q^{(t)}}(1 - {\theta ^{(t)}}) + (1 - {p^{(t)}}){{(1 - {\theta ^{(t)}})}^2}}},}&{{Z_i} = 0} \\ {\dfrac{{(1 - {p^{(t)}}){\theta ^{(t)}}{{(1 - {\theta ^{(t)}})}^2} + (1 - {p^{(t)}}){q^{(t)}}{\theta ^{(t)}}}}{{{p^{(t)}}(1 - {q^{(t)}})(1 - {\theta ^{(t)}}) + (1 - {p^{(t)}}){\theta ^{(t)}}{{(1 - {\theta ^{(t)}})}^2}}},}&{{Z_i} = 1} \\ {{Z_i},}&{{Z_i} = 2,3 \cdots } \end{array}} \right. \end{split}$

4）M步

当 $\dfrac{{\partial l(p|B,X,Y,Z)}}{{\partial p}} = 0$ ， $\dfrac{{\partial l(q|B,X,Y,Z)}}{{\partial q}} = 0$ ， $\dfrac{{\partial l(\theta |B,X,Y,Z)}}{{\partial \theta }} = 0$ 时，分别得到

$\begin{split} & {p^{(t + 1)}} = \dfrac{{\displaystyle\sum\limits_{i = 1}^n {B_i^{(t + 1)}} }}{n}\text{，}{q^{(t + 1)}} = \dfrac{{\displaystyle\sum\limits_{i = 1}^n {B_i^{(t + 1)}X_i^{(t + 1)}} }}{n}\text{，}\\ & \quad {\theta ^{(t + 1)}} = \dfrac{{\displaystyle\sum\limits_{i = 1}^n {(1 - B_i^{(t + 1)})Y_i^{(t + 1)}} }}{{\displaystyle\sum\limits_{i = 1}^n {(1 - B_i^{(t + 1)})(Y_i^{(t + 1)} + 1)} }}。 \end{split}$

根据文献[9]可得，由EM算法产生的序列 $({p^{(k)}},{q^{(k)}},{\theta ^{(k)}})$ 收敛，其极限值为 $L(p,q,\theta |B,X,Y,Z)$ 的局部最大值点。当 $({p^{(k)}},{q^{(k)}},{\theta ^{(k)}})$ 在收敛过程中小于预先指定的常数时，迭代停止，最后一次迭代值认为是EM算法下 $(p,q,\theta )$ 的极大似然估计值。

5. 贝叶斯估计

常用的先验分布为无信息先验和共轭先验，但对于ZOIGE分布，精确的共轭先验很难获得。因此，本节采用navie flat先验作为先验分布。它是一种比较简单的先验分布，参数的取值在某个范围内是均匀分布的，即对参数没有额外的先有知识。由式(3)可得， $(p,q,\theta )$ 的后验分布为

$\begin{split} & \pi (p,q,\theta |Z) \propto L(p,q,\theta |Z) = {[pq + (1 - p)(1 - \theta )]^{{M_0}}} \cdot \\ & \quad {[p(1 - q) + (1 - p)\theta (1 - \theta )]^{{M_1}}} \cdot \\ & \quad {(1 - p)^{n - {M_0} - {M_1}}}{\theta ^M}{(1 - \theta )^{n - {M_0} - {M_1}}} \end{split}$

(9)

$\begin{split} & \pi (p,q,\theta |Z) \propto \\ & \quad \displaystyle \sum\limits_j {p^{{S_{1j}}}}{q^{{S_{2j}}}}{{(1 - p)}^{{S_{3j}}}}{{(1 - q)}^{{S_{4j}}}}\cdot\\&\quad{\theta ^{M + {S_{5j}}}}{{(1 - \theta )}^{{S_{6j}}}} \end{split}$

(10)

式中： $0 < {S_{ij}} \leqslant n$ ； $i = 1,2,3,4,5,6$ 。

从式(10)可以看出，后验分布 $\pi (p,q,\theta )$ 不是常见概率分布，虽然可以直接使用马尔可夫链蒙特卡罗(Markov chain Monte Carlo, MCMC)方法来获得后验样本，但后验分布的表达式相对复杂，导致抽样的效率降低。为提高抽样效率，本节基于数据扩充下的完全似然函数和navie flat先验，将复杂的后验分布分解成简单的乘积形式，设计Gibbs抽样算法，加快抽样过程并提高效率。由式(8)可得，基于数据扩充下 $(p,q,\theta )$ 的后验分布为

$\begin{split} & \qquad \pi (p,q,\theta |B,X,Y,Z) \propto L(p,q,\theta |B,X,Y,Z) = \\ & \quad \prod\limits_{i = 1}^n {p^{{B_i}}}{{(1 - p)}^{1 - {B_i}}}{q^{{B_i}{X_i}}}{{(1 - q)}^{{B_i}(1 - {X_i})}}{\theta ^{(1 - {B_i}){Y_i}}}\cdot\\ & \quad (1 -\theta )^{(1 - {B_i})} = {p^{\sum\limits_{i = 1}^n {{B_i}} }}(1 - p)^{n - \sum\limits_{i = 1}^n {{B_i}} }{q^{\sum\limits_{i = 1}^n {{B_i}{X_i}} }}\cdot \\ & \quad(1 - p)^{\sum\limits_{i = 1}^n {{B_i}(1 - {X_i})} }{\theta ^{\sum\limits_{i = 1}^n {(1 - {B_i}){Y_i}} }}{(1 - \theta )^{n - \sum\limits_{i = 1}^n {{B_i}} }} \end{split}$

(11)

根据式(11)可得 $(p,q,\theta )$ 的满条件概率分布。

由 $\pi (p|q,\theta ,B,X,Y,Z) \propto {p^{\sum\limits_{i = 1}^n {{B_i}} }}{(1 - p)^{n - \sum\limits_{i = 1}^n {{B_i}} }}$ ， $p$ 的满条件概率分布为

$\begin{split} & p|q,\theta ,B,X,Y,Z \sim {\rm{Beta}}\\ & \quad \left(1 +\displaystyle \sum\limits_{i = 1}^n {{B_i}} ,n + 1 -\displaystyle \sum\limits_{i = 1}^n {{B_i}} \right) \end{split}$

(12)

由 $\pi (q|p,\theta ,B,X,Y,Z) \propto {q^{\sum\limits_{i = 1}^n {{B_i}{X_i}} }}{(1 - p)^{\sum\limits_{i = 1}^n {{B_i}(1 - {X_i})} }}$ ， $q$ 的满条件概率分布为

$\begin{split} & q|p,\theta ,B,X,Y,Z \sim {\rm{Beta}}\\ & \quad \left(1 +\displaystyle \sum\limits_{i = 1}^n {{B_i}{X_i}} ,1 +\displaystyle \sum\limits_{i = 1}^n {{B_i}(1 - {X_i})} \right) \end{split}$

(13)

由 $\pi (\theta |p,q,\theta ,B,X,Y,Z) \propto {\theta ^{\sum\limits_{i = 1}^n {(1 - {B_i}){Y_i}} }}{(1 - \theta )^{n - \sum\limits_{i = 1}^n {{B_i}} }}$ ， $\theta$ 的满条件概率分布为

$\begin{split} & \theta |p,q,\theta ,B,X,Y,Z \sim {\rm{Beta}}\\ & \quad \left(1 +\displaystyle \sum\limits_{i = 1}^n {(1 - {B_i})} {Y_i},n + 1 -\displaystyle \sum\limits_{i = 1}^n {{B_i}} \right) \end{split}$

(14)

结合3.2节，设计Gibbs抽样机制，具体步骤如下。

第一，设置参数初始值 $({p^{(0)}},{q^{(0)}},{\theta ^{(0)}})$ ；

第二，对 $t = 1,2 \cdots$ ，进行迭代更新：

1）利用前一步的参数估计 $({p^{(t - 1)}},{q^{(t - 1)}},{\theta ^{(t - 1)}})$ 及条件概率分布 $\pi (B,X,Y|Z,p,q,\theta )$ ，得到样本 $(B_i^{(t)},X_i^{(t)},Y_i^{(t)})$ ；

2）由 ${\rm{Beta}}\left(1 +\displaystyle \sum\limits_{i = 1}^n {B_i^{(t)}} ,n + 1 -\displaystyle \sum\limits_{i = 1}^n {B_i^{(t)}} \right)$ 产生样本 ${p^{(t)}}$ ；

3）由 ${\rm{Beta}}\left(1 +\displaystyle \sum\limits_{i = 1}^n {B_i^{(t)}X_i^{(t)}} ,1 +\displaystyle \sum\limits_{i = 1}^n {B_i^{(t)}(1 - X_i^{(t)})} \right)$ 产生样本 ${q^{(t)}}$ ；

4）由 ${\rm{Beta}}\left(1 +\displaystyle \sum\limits_{i = 1}^n {(1 - B_i^{(t)})} Y_i^{(t)},n + 1 -\displaystyle \sum\limits_{i = 1}^n {B_i^{(t)}} \right)$ 产生样本 ${\theta ^{(t)}}$ 。

6. 数值模拟

通过数值模拟，对3种ZOIGE分布参数估计的方法进行评估。样本容量分别设为 $n = 200$ 和 $n = 500$ ，几何分布中成功的概率参数 $\theta$ 分别设为 $0.3$ 和 $0.8$ ， $p$ 的值设为 $0.3$ 和 $0.7$ ， $q$ 的值设为 $0.4$ 和 $0.6$ ，所有模拟都重复2 000次。模拟过程中，计算参数估计量的均值和均方误差见表1和表2。

表 1 参数估计量的均值

Table 1. Mean of parameter estimators

参数				MLE			MLE(EM)			Bayes
$\theta$	$p$	$q$	$n$	${\theta _M}$	${p_M}$	${q_M}$	${\theta _E}$	${p_E}$	${q_E}$	${\theta _B}$	${p_B}$	${q_B}$
0.3	0.3	0.4	200	0.306	0.278	0.386	0.324	0.338	0.415	0.275	0.274	0.364
		0.4	500	0.295	0.286	0.390	0.313	0.327	0.394	0.291	0.283	0.387
		0.6	200	0.302	0.271	0.587	0.319	0.286	0.573	0.283	0.285	0.575
		0.6	500	0.296	0.282	0.593	0.305	0.291	0.586	0.288	0.294	0.588
	0.7	0.4	200	0.282	0.675	0.378	0.323	0.662	0.361	0.281	0.656	0.365
		0.4	500	0.294	0.712	0.417	0.309	0.681	0.386	0.287	0.685	0.382
		0.6	200	0.317	0.734	0.584	0.285	0.665	0.576	0.313	0.721	0.581
		0.6	500	0.293	0.717	0.590	0.291	0.684	0.588	0.306	0.709	0.585
0.8	0.3	0.4	200	0.754	0.265	0.353	0.726	0.271	0.354	0.763	0.287	0.463
		0.4	500	0.786	0.286	0.382	0.762	0.284	0.386	0.787	0.291	0.486
		0.6	200	0.722	0.275	0.577	0.779	0.279	0.548	0.761	0.283	0.576
		0.6	500	0.781	0.282	0.582	0.786	0.282	0.588	0.792	0.292	0.587
	0.7	0.4	200	0.786	0.637	0.368	0.756	0.718	0.367	0.783	0.687	0.376
		0.4	500	0.792	0.668	0.382	0.788	0.711	0.375	0.790	0.693	0.384
		0.6	200	0.776	0.686	0.564	0.772	0.728	0.578	0.784	0.685	0.564
		0.6	500	0.783	0.695	0.589	0.792	0.708	0.590	0.791	0.694	0.593

下载: 导出CSV

| 显示表格

表 2 参数估计量的均方误差

Table 2. MSE of parameter estimators

参数				MLE			MLE(EM)			Bayes
$\theta$	$p$	$q$	$n$	${\theta _{\rm{M}}}$	${p_{\rm{M}}}$	${q_{\rm{M}}}$	${\theta _{\rm{E}}}$	${p_{\rm{E}}}$	${q_{\rm{E}}}$	${\theta _{\rm{B}}}$	${p_{\rm{B}}}$	${q_{\rm{B}}}$
0.3	0.3	0.4	200	0.124	0.275	0.299	0.082	0.123	0.093	0.072	0.142	0.118
		0.4	500	0.067	0.229	0.278	0.033	0.108	0.092	0.032	0.105	0.093
		0.6	200	0.121	0.271	0.392	0.074	0.109	0.127	0.064	0.114	0.094
		0.6	500	0.068	0.232	0.361	0.034	0.113	0.159	0.023	0.076	0.059
	0.7	0.4	200	0.212	0.313	0.198	0.136	0.112	0.126	0.060	0.244	0.064
		0.4	500	0.107	0.297	0.164	0.085	0.109	0.115	0.046	0.115	0.056
		0.6	200	0.205	0.307	0.278	0.122	0.218	0.189	0.053	0.098	0.091
		0.6	500	0.114	0.293	0.253	0.081	0.196	0.109	0.048	0.073	0.084
0.8	0.3	0.4	200	0.114	0.085	0.143	0.074	0.071	0.154	0.113	0.187	0.133
		0.4	500	0.093	0.044	0.134	0.068	0.056	0.134	0.086	0.143	0.126
		0.6	200	0.122	0.055	0.141	0.058	0.021	0.148	0.121	0.124	0.131
		0.6	500	0.091	0.038	0.122	0.023	0.018	0.129	0.091	0.117	0.127
	0.7	0.4	200	0.084	0.067	0.142	0.154	0.082	0.043	0.087	0.163	0.141
		0.4	500	0.069	0.042	0.138	0.108	0.068	0.035	0.067	0.144	0.137
		0.6	200	0.054	0.064	0.094	0.128	0.032	0.078	0.076	0.118	0.162
		0.6	500	0.049	0.034	0.089	0.113	0.029	0.062	0.057	0.115	0.149

下载: 导出CSV

| 显示表格

由表1可知，当 $\theta$ 较小时，3种估计方法在估计值的精确性方面大致相当，随着样本容量的增加，估计值越来越接近真实值；当 $\theta$ 较大时，极大似然估计与EM算法下的极大似然估计表现较好，而贝叶斯估计方法表现欠佳。当 $p$ 较小时，3种方法下 $\theta$ 的估计值在精确性方面要高于 $p$ 较大时的情况，这是因为当 $p$ 较大时，样本来源于几何分布的概率较小，影响了估计效果。

由表2可知，当 $\theta$ 较小时，贝叶斯方法比其他两种估计方法更稳定；当 $\theta$ 较大时，其他两种方法稳定性相当，而贝叶斯方法则没有其余两种估计方法的稳定性好。

7. 实例分析

本节对R包NMMAPS中1994年底特律每日交通事故导致死亡的数据集进行研究，该数据集共有365个观测数据点，拟合结果如图1所示，0和1的观测频数较高，表现出明显的0−1膨胀现象。然而传统的泊松分布和负二项分布没有很好的拟合效果，故采用本文模型进行拟合。在拟合过程中，对于EM算法下的极大似然估计，E步和M步的迭代次数为500次，在贝叶斯方法中，隐变量的更新次数为3 000次，burn-in为2 000次。参数的点估计见表3。参数的拟合值见表4。可以看出，从比较AIC及拟合值与观测值的接近程度来看，贝叶斯方法的估计相较于其他两种方法效果更佳。

图 1 1994年底特律交通事故死亡数量的拟合频数

Figure 1. Fitted frequency of traffic accident deaths in Detroit in 1994

下载: 全尺寸图片幻灯片

表 3 3种方法下参数的点估计

Table 3. Point estimation of parameters under three methods

参数	MLE	MLE(EM)	Bayes
$\theta$	0.0131	0.0131	0.0127
$p$	0.6830	0.6804	0.6616
$q$	0.6485	0.6490	0.6515
AIC	670.0665	668.0072	660.273

下载: 导出CSV

| 显示表格

表 4 3种方法下观测值与拟合值的比较

Table 4. Comparison of observed and fitted values under three methods

观测值	观测频数	MLE	MLE(EM)	Bayes
0	181	180	180	181
1	122	120	120	121
2	28	31	31	29
3	25	27	27	28
4	5	5	6	4
5	2	2	1	1
6	1	0	0	1
7	1	0	0	0

下载: 导出CSV

| 显示表格

8. 结语

鉴于样本数据中存在较多的零、一观测值，本研究建立了0−1膨胀几何分布模型，采用极大似然估计法、EM算法下的极大似然估计法和贝叶斯方法进行参数估计。为解释样本数据中零一膨胀现象，引入隐变量，将复杂的似然函数转化为简单的乘积形式，便于进行理论分析与试验模拟。在参数估计方面，研究表明，EM算法下的极大似然估计和贝叶斯估计表现更好，随着样本量的增加，得到的参数估计更加接近真实值，估计结果更为无偏和有效。在未来的研究中，可继续探索其他的研究方法，进一步提高对膨胀数据的建模和预测能力。

图 1 1994年底特律交通事故死亡数量的拟合频数

Figure 1. Fitted frequency of traffic accident deaths in Detroit in 1994

下载: 全尺寸图片幻灯片

表 1 参数估计量的均值

Table 1. Mean of parameter estimators

参数				MLE			MLE(EM)			Bayes
$\theta$	$p$	$q$	$n$	${\theta _M}$	${p_M}$	${q_M}$	${\theta _E}$	${p_E}$	${q_E}$	${\theta _B}$	${p_B}$	${q_B}$
0.3	0.3	0.4	200	0.306	0.278	0.386	0.324	0.338	0.415	0.275	0.274	0.364
		0.4	500	0.295	0.286	0.390	0.313	0.327	0.394	0.291	0.283	0.387
		0.6	200	0.302	0.271	0.587	0.319	0.286	0.573	0.283	0.285	0.575
		0.6	500	0.296	0.282	0.593	0.305	0.291	0.586	0.288	0.294	0.588
	0.7	0.4	200	0.282	0.675	0.378	0.323	0.662	0.361	0.281	0.656	0.365
		0.4	500	0.294	0.712	0.417	0.309	0.681	0.386	0.287	0.685	0.382
		0.6	200	0.317	0.734	0.584	0.285	0.665	0.576	0.313	0.721	0.581
		0.6	500	0.293	0.717	0.590	0.291	0.684	0.588	0.306	0.709	0.585
0.8	0.3	0.4	200	0.754	0.265	0.353	0.726	0.271	0.354	0.763	0.287	0.463
		0.4	500	0.786	0.286	0.382	0.762	0.284	0.386	0.787	0.291	0.486
		0.6	200	0.722	0.275	0.577	0.779	0.279	0.548	0.761	0.283	0.576
		0.6	500	0.781	0.282	0.582	0.786	0.282	0.588	0.792	0.292	0.587
	0.7	0.4	200	0.786	0.637	0.368	0.756	0.718	0.367	0.783	0.687	0.376
		0.4	500	0.792	0.668	0.382	0.788	0.711	0.375	0.790	0.693	0.384
		0.6	200	0.776	0.686	0.564	0.772	0.728	0.578	0.784	0.685	0.564
		0.6	500	0.783	0.695	0.589	0.792	0.708	0.590	0.791	0.694	0.593

下载: 导出CSV

表 2 参数估计量的均方误差

Table 2. MSE of parameter estimators

参数				MLE			MLE(EM)			Bayes
$\theta$	$p$	$q$	$n$	${\theta _{\rm{M}}}$	${p_{\rm{M}}}$	${q_{\rm{M}}}$	${\theta _{\rm{E}}}$	${p_{\rm{E}}}$	${q_{\rm{E}}}$	${\theta _{\rm{B}}}$	${p_{\rm{B}}}$	${q_{\rm{B}}}$
0.3	0.3	0.4	200	0.124	0.275	0.299	0.082	0.123	0.093	0.072	0.142	0.118
		0.4	500	0.067	0.229	0.278	0.033	0.108	0.092	0.032	0.105	0.093
		0.6	200	0.121	0.271	0.392	0.074	0.109	0.127	0.064	0.114	0.094
		0.6	500	0.068	0.232	0.361	0.034	0.113	0.159	0.023	0.076	0.059
	0.7	0.4	200	0.212	0.313	0.198	0.136	0.112	0.126	0.060	0.244	0.064
		0.4	500	0.107	0.297	0.164	0.085	0.109	0.115	0.046	0.115	0.056
		0.6	200	0.205	0.307	0.278	0.122	0.218	0.189	0.053	0.098	0.091
		0.6	500	0.114	0.293	0.253	0.081	0.196	0.109	0.048	0.073	0.084
0.8	0.3	0.4	200	0.114	0.085	0.143	0.074	0.071	0.154	0.113	0.187	0.133
		0.4	500	0.093	0.044	0.134	0.068	0.056	0.134	0.086	0.143	0.126
		0.6	200	0.122	0.055	0.141	0.058	0.021	0.148	0.121	0.124	0.131
		0.6	500	0.091	0.038	0.122	0.023	0.018	0.129	0.091	0.117	0.127
	0.7	0.4	200	0.084	0.067	0.142	0.154	0.082	0.043	0.087	0.163	0.141
		0.4	500	0.069	0.042	0.138	0.108	0.068	0.035	0.067	0.144	0.137
		0.6	200	0.054	0.064	0.094	0.128	0.032	0.078	0.076	0.118	0.162
		0.6	500	0.049	0.034	0.089	0.113	0.029	0.062	0.057	0.115	0.149

下载: 导出CSV

表 3 3种方法下参数的点估计

Table 3. Point estimation of parameters under three methods

参数	MLE	MLE(EM)	Bayes
$\theta$	0.0131	0.0131	0.0127
$p$	0.6830	0.6804	0.6616
$q$	0.6485	0.6490	0.6515
AIC	670.0665	668.0072	660.273

下载: 导出CSV

表 4 3种方法下观测值与拟合值的比较

Table 4. Comparison of observed and fitted values under three methods

观测值	观测频数	MLE	MLE(EM)	Bayes
0	181	180	180	181
1	122	120	120	121
2	28	31	31	29
3	25	27	27	28
4	5	5	6	4
5	2	2	1	1
6	1	0	0	1
7	1	0	0	0

下载: 导出CSV

[1]	张良超, 周金亮, 温利民. 零膨胀泊松模型中风险参数的贝叶斯估计[J] . 江西师范大学学报(自然科学版),2020,44(3):269 − 274. doi: 10.16357/j.cnki.issn1000-5862.2020.03.10
[2]	肖翔, 刘福窑. 零膨胀几何分布的参数估计[J] . 上海工程技术大学学报,2018,32(3):267 − 271,277. doi: 10.3969/j.issn.1009-444X.2018.03.013
[3]	ZHANG C, TIAN G L, NG K W. Properties of the zero-and-one inflated poisson distribution and likelihood-based inference methods[J] . Statistics and Its Interface,2016,9(1):11 − 32. doi: 10.4310/SII.2016.v9.n1.a2
[4]	TANG Y C, LIU W C, XU A C. Statistical inference for zero-and-one-inflated poisson models[J] . Statistical Theory and Related Fields,2017,1(2):216 − 226. doi: 10.1080/24754269.2017.1400419
[5]	JORNSATIAN C, BODHISUWAN W. Zero-one inflated negative binomial-beta exponential distribution for count data with many zeros and ones[J] . Communications in Statistics: Theory and Methods,2022,51(24):8517 − 8531. doi: 10.1080/03610926.2021.1898642
[6]	TAJUDDIN R R, ISMAIL N, IBRAHIM K, et al. A new zero–one-inflated poisson–lindley distribution for modelling overdispersed count data[J] . Bulletin of the Malaysian Mathematical Sciences Society,2022,45(1):21 − 35.
[7]	XIAO X, TANG Y C, XU A C, et al. Bayesian inference for zero-and-one-inflated geometric distribution regression model using Pólya-Gamma latent variables[J] . Communications in Statistics - Theory and Methods,2020,49(15):3730 − 3743. doi: 10.1080/03610926.2019.1709647
[8]	肖翔, 古晞. 0–1膨胀几何分布的客观贝叶斯分析[J] . 上海工程技术大学学报,2021,35(3):266 − 271.
[9]	WU C F J. On the convergence properties of the EM algorithm[J] . The Annals of Statistics,1983,11(1):95 − 103.

姓名
邮箱
手机号码
标题
留言内容
验证码

留言板