混合幂级数分布及其应用

Mixture Power Series Distribution and Its Application

摘要: 幂级数分布是一类广泛的概率分布，包括了常用的二项分布、Poisson分布等多种离散型概率分布。由若干个幂级数分布构成的混合幂级数分布，能够灵活地拟合如零膨胀、0~k膨胀以及其他多种膨胀类型的计数数据。本文基于EM算法给出了混合幂级数分布的参数估计，并对混合幂级数分布的特例0~1膨胀Poisson分布进行了说明，通过实例分析，说明混合幂级数分布及估计方法是适用的。

Abstract: Power series distribution is a kind of extensive probability distribution, including the commonly used binomial distribution, Poisson distribution and other many discrete distributions. The mixed power series distribution, which is composed of multiple power series distributions, can flexibly fit the counting data of zero inflated, 0~k inflated and other many of inflated types. In this paper, the parameter estimation of mixed power series distribution is given based on EM algorithm, and the special case of mixed power series distribution 0~1 inflated Poisson distribution is explained. Through case analysis, the mixed power series distribution and estimation method are proved to be applicable.

1. 引言

幂级数分布最初是由Noak [1] 在研究了一类具有共同特性的离散型随机变量分布的基础上，给出了幂级数分布的定义，同时也研究了关于幂级数分布矩的性质，并对幂级数分布中的一些重要特例，如二项分布、Poisson分布、负二项分布等分别进行了说明。幂级数分布属于离散型概率分布族，目前已推广到了广义幂级数分布 [2] ，如广义Poisson分布和广义负二项分布等都属于这一类型的分布。幂级数分布适用于对计数数据的统计分析，而计数数据是一类重要的数据类型，广泛存在于医学、保险精算、人口学等多个研究领域中。

在实际对计数数据的统计过程中，经常会出现观测值在某些特定数值上频数过高的现象，例如在某段时间内汽车保险索赔次数、人工流产的次数等等，就是数据中0出现的频数明显过高，高于预期分布中0的期望值，即所谓的零膨胀数据。对于该类型的数据统计方法的研究，解锋昌、韦博成和林金官 [3] 有比较系统的论述。

除了零膨胀现象外，还会遇到某个特定的正整数值k大量出现的情况，如某个地区居民一年内看牙医的次数，观测数据中大量出现了0和1，这可能是因为许多人没有养成定期看牙医的良好习惯，或是牙齿有问题去看牙医之后，短期之内就不再去看牙医。为了更好地拟合此类数据，Lin和Tsai [4] 提出了0~k膨胀的模型。

在目前已有的研究中，对零膨胀数据或0~k膨胀数据进行分析时，模型大多选用的是Poisson分布或负二项分布，是否存在其他分布更适合这种数据，是一个值得探讨的问题。

2. 混合幂级数分布

2.1. 幂级数分布

由幂级数展开式 $g\left(\theta \right)={\displaystyle \underset{l=0}{\overset{\infty }{\sum }}a\left(l\right){\theta }^l}$ ，其中， $a\left(l\right)$ 为常数序列，从而有

${\displaystyle \underset{l=0}{\overset{\infty }{\sum }}\frac{a\left(l\right){\theta }^l}{g\left(\theta \right)}}=1$

故定义离散型随机变量X的概率分布律为

$P\left(X=l\right)=\frac{a\left(l\right){\theta }^l}{g(\theta )}$ ， $l=0,1,2，\cdots$ 且 $\theta >0$ ， $a_l>0$ ，

称X为幂级数分布。

如在幂级数分布定义中， $a_l=\frac{1}{l!}$ ， $g\left(\theta \right)={\text{e}}^{\theta }$ ，显然X为Possion分布，记为 $X~\text{Possion}\left(\theta \right)$ 。又如取 $a_l=C_n^l$ ， $g\left(\theta \right)={\left(1+\theta \right)}^n={\displaystyle \underset{l=0}{\overset{\infty }{\sum }}{C}_n^l}{\theta }^l$ 时，X的概率分布律为

$P\left(X=l\right)=C_n^l\frac{{\theta }^l}{{\left(1+\theta \right)}^n}=C_n^l{\left(\frac{\theta }{1+\theta }\right)}^l{\left(\frac{1}{1+\theta }\right)}^{n-l}$ ， $l=0,1,2，\cdots ,n$ ，

是二项分布，记为 $X~\text{b}\left(n,\frac{\theta }{1+\theta }\right)$ 。

2.2. 混合幂级数分布

设离散型随机变量Y，概率分布律为

$f\left(y|\theta \right)={\displaystyle \underset{k=1}{\overset{m}{\sum }}{\alpha }_k}{f}_k\left(y|{\theta }_k\right)={\alpha }_1{f}_1\left(y|{\theta }_1\right)+{\alpha }_2{f}_2\left(y|{\theta }_2\right)+\cdots +{\alpha }_m{f}_m\left(y|{\theta }_m\right)$ (1)

其中， $f_k\left(y|{\theta }_k\right)=\frac{a_k\left(y\right){\theta }_k^y}{g_k\left({\theta }_k\right)}$ ， $\theta =\left({\theta }_1,\cdots ,{\theta }_m\right)$ ， $y=0,1,\cdots$ ， ${\alpha }_k$ 为系数且 ${\alpha }_k>0$ ， ${\displaystyle \underset{k=1}{\overset{m}{\sum }}{\alpha }_k}=1$ ，称为混合幂级数分布。

混合幂级数分布的概率母函数为

$G_Y\left(s\right)=E\left(s^Y\right)={\displaystyle \underset{k=1}{\overset{m}{\sum }}{\alpha }_k{\displaystyle \underset{y=0}{\overset{\infty }{\sum }}{f}_k\left(y|{\theta }_k\right)s^y}}={\displaystyle \underset{k=1}{\overset{m}{\sum }}{\alpha }_k{\displaystyle \underset{l=0}{\overset{\infty }{\sum }}\frac{a_k\left(l\right){\left({\theta }_{k}s\right)}^l}{g_k(\; \theta \; k\; )}}}$

通过概率母函数可以求得混合幂级数分布的各阶矩。

3. 混合幂级数分布的参数估计

若构成混合分布的幂级数分布较多，则需要估计的参数也多，用矩估计法是比较困难的。因此，对混合幂级数分布参数估计主要讨论EM算法 [5] 。

设观测 $Y_1,Y_2,\cdots ,Y_n$ 是来自总体分布如式(1)的混合幂级数分布的样本，引入隐变量 $Z_{ki}$ ，当 $Y_i$ 来自分布 $f_k\left(y|{\theta }_k\right)$ 时，有 $f\left(Z_{ki}=1\right)={\alpha }_k$ ， $k=1,2,\cdots ,m$ ，记

$y=\left(y_1,\cdots ,y_n\right)$ ， $Z_k=\left(Z_{k1},\cdots ,Z_{kn}\right)$ ， $\vartheta =\left({\theta }_1,\cdots ,{\theta }_m,{\alpha }_1,\cdots ,{\alpha }_m\right)$ ，

完全数据的似然函数为

$L\left(\vartheta |y,z_k\right)={\displaystyle \underset{k=1}{\overset{m}{\prod }}{\displaystyle \underset{i=1}{\overset{n}{\prod }}{\left[{\alpha }_k{f}_k\left(y_i|{\theta }_k\right)\right]}^{z_{ki}}}}$ ，

对数似然函数为

$l\left(\vartheta |y,z_k\right)=\log L\left(\vartheta |y,z_k\right)={\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{z}_{ki}}\left[\log {\alpha }_k{f}_k\left({\theta }_k|y_i\right)\right]}$

E步：

${\hat{z}}_{ki}=E\left(z_{ki}|y_i,{\theta }_k\right)=f_k\left(z_{ki}|y_i,{\theta }_k\right)=\frac{f_k\left(z_{ki},y_i|{\theta }_k\right)}{{\displaystyle \underset{k=1}{\overset{m}{\sum }}{f}_k\left(z_{ki},y_i|{\theta }_k\right)}}=\frac{{\alpha }_k{f}_k\left(y_i|{\theta }_k\right)}{{\displaystyle \underset{k=1}{\overset{m}{\sum }}{\alpha }_k{f}_k\left(y_i|{\theta }_k\right)}}=\frac{{\alpha }_k\frac{a_k\left(y_i\right){\theta }_k^{y_{{}_i}}}{g_k\left({\theta }_k\right)}}{{\displaystyle \underset{k=1}{\overset{m}{\sum }}{\alpha }_k\frac{a_k\left(y_i\right){\theta }_k^{y_{{}_i}}}{g_k(\; \theta \; k\; )}}}$

令

$Q\left(\vartheta ,{\vartheta }^{\left(s\right)}\right)=E\left\{{\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{z}_{ki}}\left[\log {\alpha }_k{f}_k\left({\theta }_k|y_i\right)\right]}\right\}$

其中， ${\vartheta }^{\left(s\right)}$ 表示第s次迭代的各参数，则

$\begin{array}{c}Q\left(\vartheta ,{\vartheta }^{\left(s\right)}\right)={\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\hat{z}}_{ki}}\left[\log {\alpha }_k+\log f_k\left(y_i|{\theta }_k\right)\right]}\\ ={\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\hat{z}}_{ki}}\left[\log {\alpha }_k+\log a_k\left(y_i\right)+y_i\log {\theta }_k-\log g_k\left({\theta }_k\right)\right]}\end{array}$

M步：求函数 $Q\left(\vartheta ,{\vartheta }^{\left(s\right)}\right)$ 的极大值点，计算的迭代公式为

$\{\begin{array}{l}{\hat{\alpha }}_k^{\left(s+1\right)}=\frac{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\hat{z}}^{\left(s\right)}{}_{jk}}}{{\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\hat{z}}^{\left(s\right)}{}_{jk}}}}\\ {\hat{\theta }}_k^{\left(s+1\right)}=\frac{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\hat{z}}^{\left(s\right)}{}_{jk}{y}_j{f}_k\left({\theta }_k^{\left(s\right)}\right)}}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\hat{z}}^{\left(s\right)}{}_{jk}{f}_k\left({\theta }_k^{\left(s\right)}\right)}}\end{array}$ (2)

重复上述步骤，直到两次迭代值之差小于给定值为止。通过幂级数分布的一些特例，如二项分布、Poisson分布、几何分布、负二项分布等可以构造出各种形式灵活的混合幂级数分布。

4. 混合幂级数分布的特例及应用

本节主要通过0~1分布和Poisson分布构造出混合幂级数分布的特例0~1膨胀Poisson分布。

4.1. 0~1膨胀Poisson分布

设随机变量 $X_1~\text{b}\left(1,\frac{{\theta }_1}{1+{\theta }_1}\right)$ ， $X_2~\text{Possion}\left({\theta }_2\right)$ ， $Z=\{\begin{array}{l}1当Y=X_1时\\ 0当Y=X_2时\end{array}$ ，

且 $f\left(Z=1\right)=\alpha$ ， $f\left(Z=0\right)=1-\alpha$ ，即 $Y=Z{X}_1+\left(1-Z\right)X_2$ ，Y的概率分布律为

$P\left(Y=k\right)=\{\begin{array}{l}\alpha \frac{1}{1+{\theta }_1}+\left(1-\alpha \right){\text{e}}^{-{\theta }_2}k=0\\ \alpha \frac{{\theta }_1}{1+{\theta }_1}+\left(1-\alpha \right){\theta }_2{\text{e}}^{-{\theta }_2}k=1\\ \left(1-\alpha \right)\frac{{\theta }_2^k}{k!}{\text{e}}^{-{\theta }_2}k=2,3,\cdots \end{array}$ (3)

称Y服从0~1膨胀Poisson分布。

4.2. 0~1膨胀Poisson分布的参数估计

4.2.1. 矩估计

0~1膨胀Poisson分布的概率母函数为

$G_Y\left(t\right)=E\left(t^Y\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{t}^k}p\left(Y=k\right)=\alpha \left(\frac{1}{1+{\theta }_1}+\frac{t{\theta }_1}{1+{\theta }_1}\right)+\left(1-\alpha \right){\text{e}}^{{\theta }_2\left(t-1\right)}$ ，

求它的1阶、2阶、3阶导数分别为

${G^{\prime }}_Y\left(1\right)=\alpha \frac{{\theta }_1}{1+{\theta }_1}+\left(1-\alpha \right){\theta }_2$ ， ${G^{″}}_Y\left(1\right)=\left(1-\alpha \right){\theta }_2^2$ ， ${G^{‴}}_Y\left(1\right)=\left(1-\alpha \right){\theta }_2^3$ ，

因此

$E\left(Y\right)=G^{\prime }\left(1\right)=\alpha \frac{{\theta }_1}{1+{\theta }_1}+\left(1-\alpha \right){\theta }_2$ ，

$E\left(Y^2\right)=G^{″}\left(1\right)+G^{\prime }\left(1\right)=\left(1-\alpha \right){\theta }_2^2+\left(1-\alpha \right){\theta }_2+\alpha \frac{{\theta }_1}{1+{\theta }_1}$ ，

$E\left(Y^3\right)=G^{‴}\left(1\right)+3G^{″}\left(1\right)+G^{\prime }\left(1\right)=\left(1-\alpha \right){\theta }_2^3+3\left(1-\alpha \right){\theta }_2^2+\left(1-\alpha \right){\theta }_2+\alpha \frac{{\theta }_1}{1+{\theta }_1}$ 。

设 $Y_1,Y_2,\cdots ,Y_n$ 是从式(3)的总体Y中抽取的样本，样本的1阶、2阶和3阶原点矩分别记为

$m_1=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{Y}_i}$ ， $m_2=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{Y}_i^2}$ ， $m_3=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{Y}_i^3}$ ，

由矩估计的原理，可得

${\hat{\theta }}_2=\frac{m_3-3m_2+2m_1}{m_2-m_1}$ ， $\hat{\alpha }=1-\frac{{\left(m_2-m_1\right)}^3}{{\left(m_3-3m_2+2m_1\right)}^2}$ (4)

进而得

${\hat{\theta }}_1=\frac{m_1-\left(1-\hat{\alpha }\right){\hat{\theta }}_2}{\hat{\alpha }-m_1+\left(1-\hat{\alpha }\right){\hat{\theta }}_2}$ (5)

4.2.2. EM估计

总体Y为0~1膨胀Poisson分布，在混合幂级数分布式(1)的定义，取 $m=2$ ，只需引入隐变量为 $Z_i$ ，当 $Y_i$ 来自分布 $Y_i~\text{b}\left(1,\frac{{\theta }_1}{1+{\theta }_1}\right)$ 时，有 $f\left(Z_i=1\right)=\alpha$ ，当 $Y_i$ 来自分布 $\text{Possion}\left({\theta }_2\right)$ 时，有 $f\left(Z_i=0\right)=1-\alpha$ ， $i=1,2,\cdots ,n$ 。

记 $z=\left(z_1,\cdots ,z_n\right)$ ， $y=\left(y_1,\cdots ,y_n\right)$ ，完全数据的似然函数为

$L\left(\theta |y,z\right)={\displaystyle \underset{i=1}{\overset{n}{\prod }}{\left[\alpha \frac{{\theta }_1^{y_i}}{1+{\theta }_1}\right]}^{z_i}}{\left[\left(1-\alpha \right)\frac{{\theta }_2^{y_i}}{x_i!}{\text{e}}^{-{\theta }_2}\right]}^{1-z_i}$ ，

对数似然函数为

$\mathcal{l}\left(\theta |y,z\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left\{z_i\left[\ln \alpha +y_i\ln {\theta }_1-\ln \left(1+{\theta }_1\right)\right]+\left(1-z_i\right)\left[\ln \left(1-\alpha \right)+y_i\ln {\theta }_2-{\theta }_2-\ln y_i!\right]\right\}}$ ，

然后根据本文第2部分EM算法中的迭代公式(2)可以得到参数估计结果。

事实上，对于0~1膨胀Poisson分布，在E步中得到的 ${\hat{z}}_i$ 更简单，只有三种情形。具体地

当 $y_i=0$ 时， ${\hat{z}}_i=\frac{\frac{\alpha }{1+{\theta }_1}}{\frac{\alpha }{1+{\theta }_1}+\left(1-\alpha \right){\text{e}}^{-{\theta }_2}}=\frac{1}{1+\left(1+{\theta }_1\right)\left(1/\alpha -1\right){\text{e}}^{-{\theta }_2}}$ ，记为 ${\tilde{p}}_0$ ；

当 $y_i=1$ 时， ${\hat{z}}_i=\frac{\frac{\alpha {\theta }_1}{1+{\theta }_1}}{\frac{\alpha {\theta }_1}{1+{\theta }_1}+\left(1-\alpha \right){\theta }_2{\text{e}}^{-{\theta }_2}}=\frac{1}{1+\left(1+1/{\theta }_1\right)\left(1/\alpha -1\right){\theta }_2{\text{e}}^{-{\theta }_2}}$ ，记为 ${\tilde{p}}_1$ ；

当 $y_i=2,3$ 时， ${\hat{z}}_i=0$ 。

记 $n_0$ 为样本中 $y_i=0$ 的个数， $n_1$ 为样本中 $y_i=1$ 的个数， ${\theta }^{\left(0\right)}$ 为初始参数值，故在M步中，

$\begin{array}{c}Q\left(\theta |y,{\theta }^{\left(0\right)}\right)=\left\{{\tilde{p}}_0\left[\ln \alpha -\ln \left(1+{\theta }_1\right)\right]+\left(1-{\tilde{p}}_0\right)\left[\ln \left(1-\alpha \right)-{\theta }_2\right]\right\}{n}_0\\ +\left\{{\tilde{p}}_1\left(\ln \alpha +\ln \frac{{\theta }_1}{1+{\theta }_1}\right)+\left(1-{\tilde{p}}_1\right)\left[\ln \left(1-\alpha \right)+\ln {\theta }_2-{\theta }_2\right]\right\}{n}_1\\ +{\displaystyle \underset{y_i>2}{\sum }\left[\ln \left(1-\alpha \right)+y_i\ln {\theta }_2-\ln y_i-{\theta }_2\right]}\end{array}$

从而，得

$\{\begin{array}{l}{\hat{\alpha }}^{\left(s+1\right)}=\frac{n_0{\tilde{p}}_0^{\left(s\right)}+n_1{\tilde{p}}_1^{\left(s\right)}}{n}\\ {\hat{\theta }}_1^{\left(s+1\right)}=\frac{n_1{\tilde{p}}_1^{\left(s\right)}}{n_0{\tilde{p}}_0^{\left(s\right)}}\\ {\hat{\theta }}_2^{\left(s+1\right)}=\frac{{\displaystyle \sum y_i-n_1{\tilde{p}}_1^{\left(s\right)}}}{n-n_0{\tilde{p}}_0^{\left(s\right)}-n_1{\tilde{p}}_1^{\left(s\right)}}\end{array}$ (6)

其中， ${\tilde{p}}_0^{\left(s\right)}$ 、 ${\tilde{p}}_1^{\left(s\right)}$ 为第s步迭代值， ${\hat{\alpha }}^{\left(s+1\right)}$ 、 ${\hat{\theta }}_1^{\left(s+1\right)}$ 、 ${\hat{\theta }}_2^{\left(s+1\right)}$ 为第 $s+1$ 步的迭代值。通过上述迭代，得到参数的最大似然估计值。

4.3. 实例分析

本小节分析的数据集是新西兰白兔出产死亡的数据，最初是由Morgan等 [6] 研究得分检验时使用过。表1的第1列和第2列列出了新西兰白兔每窝死产数量和观测频数。数据分别用零膨胀Poisson分布(ZIP)和0~1膨胀Poisson (ZIOP)分布进行拟合的，参数估计的方法采用矩估计(ME)和基于EM算法的最大似然估计(MLE)。

Table 1. Fitting results of stillbirth numbers in New Zealand white rabbits

表1. 新西兰白兔死产数量数据的拟合结果

ZOIP模型的参数的矩估计通过式(4)和式(5)，计算结果分别列于表1的第4列中的第12行至第14行，最大似然估计通过式(6)得到，计算结果分别列于表1的第6列中的第12行至第14行。作为对照，将ZIP模型的参数的矩估计和极大似然估计结果分别列于表1的第3列和第5列中的第12行和第14行。

根据模型参数估计值拟合的新西兰白兔死产的频数列于表1上半部分，ZIP模型和ZOIP模型的拟合优度检验统计量列于表1的下半部分。由于ZIP模型包括2个参数，ZOIP模型3个参数，拟合优度检验中ZIP模型自由度为 $7-2-1=4$ ，ZOIP模型自由度为 $7-3-1=3$ 。拟合优度检验结果表明，只有基于EM算法得到的ZOIP模型的拟合优度检验P值大于0.05，故可以认为基于EM算法的ZOIP模型拟合该数据集是适合的。

5. 讨论

事实上，为了适应更多类型的数据，需要尝试选择多种分布来拟合数据。因此，对幂级数分布及相关内容的研究具有重要意义。如在幂级数分布基础上可以选择加入一个比例参数，构成零膨胀幂级数分布，也可以选择若干个幂级数分布的构成混合幂级数分布。混合幂级数分布可以构造出很多灵活的模型，包括零膨胀或0~k膨胀以及多种膨胀类型的计数数据，拓宽了幂级数分布的应用范围，也为多点膨胀计数模型提供了一个新思路。

基金项目

本研究由2022年度辽宁省教育厅高校基本科研项目(LJKMZ20221412)和2022年度辽宁省研究生教育教学改革研究项目(2022-180-39510165)资助。

NOTES

^*通讯作者。