Lindley分布参数变点的贝叶斯估计

Bayesian Estimation of Parameter Change Points of Lindley Distribution

摘要: 利用贝叶斯方法研究了Lindley分布参数存在变点的参数估计问题，给出Lindley分布的变点模型，对参数选取无信息先验分布和伽玛分布两种情况，分别求出各参数的满条件分布，并通过R软件做随机模拟，得出各参数的MC误差都小于2%，且区间估计效果理想，表明通过贝叶斯估计研究各参数的估计值是有效的。

Abstract: The parameter estimation problem of Lindley distribution with change points is studied by using Bayesian method. The change point model of Lindley distribution is given. The full conditional distribution of each parameter is calculated for the two cases of no information prior distribution and gamma distribution when the parameters are selected. The random simulation by R software shows that the MC error of each parameter is less than 2%, and the interval estimation effect is ideal. It shows that the estimation of each parameter by Bayesian estimation is effective.

1. 引言

Lindley分布是可靠性研究中的一个重要分布，某些寿命数据可通过Lindley模型达到更好的拟合效果 [1]。Krishna等在逐步II型右删失数据下，采用极大似然方法和贝叶斯方法研究了Lindley分布的可靠性 [2]。杨冬霞等分别在完全数据、逐步I型区间删失数据，逐步II型删失数据以及定数截尾样本下研究了Lindley分布的参数估计问题 [3] [4] [5]；范梓淼等分别讨论了在NA随机样本序列和独立同分布样本下Lindley分布参数的经验贝叶斯检验函数问题 [6] [7]；龙兵分析了Lindley分布参数的区间估计和假设检验问题 [8]；近几年，变点问题也是统计方向研究的一个热点问题。何朝兵等在左截断右删失数据下对指数分布多变点模型进行了参数估计 [9]。沙雪云等利用贝叶斯方法研究了Lomax分布形状参数变点的估计模型 [10]；程静等用极大似然估计和贝叶斯估计讨论了两种分布的单变点问题 [11] [12]。关于Lindley分布参数变点问题的研究较少，本文给出了Lindley分布的变点模型，在叙述了解决多变点模型问题的具体步骤后，主要研究Lindley分布参数的单变点模型，分别在无信息先验分布和伽玛分布为先验分布的条件下，利用贝叶斯估计研究参数和变点位置，并通过R软件进行随机模拟。结果显示：各参数的估计值和真实值之间的MC误差较小，表明其估计值的效果较为理想。

2. Lindley分布变点模型

设随机变量 服从参数为 $\theta$ 的Lindley分布，则分布函数和密度函数如下：

$F\left(x;\theta \right)=1-\left(1+\frac{\theta }{\theta +1}x\right)\exp \left(-\theta x\right),x>0.$

$f\left(x;\theta \right)=\frac{{\theta }^2}{\theta +1}\left(1+x\right)\exp \left(-\theta x\right),x>0.$

其中参数 $\theta >0$。

Lindley分布多变点模型为：

$X_i~\{\begin{array}{l}\text{Lindley}\left({\theta }_1\right),i=1,2,\cdots ,k_1\\ \text{Lindley}\left({\theta }_2\right),i=k_1+1,\cdots ,k_2\\ ⋮\\ \text{Lindley}\left({\theta }_{m+1}\right),i=k_m+1,\cdots ,n.\end{array}$

其中 ${\theta }_1,{\theta }_2,\cdots ,{\theta }_{m+1}$ 两两不等，m是变点个数， $k_1,k_2,\cdots ,k_m$ (满足 $1\le k_1<k_2<\cdots <k_m\le n-1$ )是需要估计的变点位置。通过二分分段法来解决多变点的问题的具体步骤为：先确定Lindley分布的序列S中是否存在单变点，如果没有，则序列S中无变点；如果存在单变点，此变点将Lindley分布的序列S拆分成两个子序列，再次确定两个子序列中是否存在单变点，重复上述步骤，直至所有子序列中识别不到变点为止。

设随机变量 $X_i\left(i=1,2,\cdots ,n\right)$ 相互独立且满足

$X_i~\{\begin{array}{l}\text{Lindley}\left({\theta }_1\right),i=1,2,\cdots ,k,\\ \text{Lindley}\left({\theta }_2\right),i=k+1,\cdots ,n.\end{array}$

其中参数 ${\theta }_1,{\theta }_2>0,1\le k\le n-1$ 且 $k\in N^{+}$， $k,{\theta }_1,{\theta }_2$ 均未知，当 ${\theta }_1\ne {\theta }_2$ 时k就是要讨论的变点，此模型只含有一个变点，称其为Lindley分布的单变点模型。

下文确定各参数的贝叶斯估计，对 取无信息先验分布： $\pi \left(k\right)=\frac{1}{n-1}$，对参数 ${\theta }_1,{\theta }_2$ 分别取无信息先验分布和伽玛分布后，再对变点k和参数 ${\theta }_1,{\theta }_2$ 做贝叶斯估计。

3. Lindley分布参数的贝叶斯估计

当 ${\theta }_1\ne {\theta }_2$ 时，设k是变点，故此变点问题的似然函数为

$\begin{array}{c}L\left(k,{\theta }_1,{\theta }_2|x\right)=\frac{{\theta }_1^2}{{\theta }_1+1}\left(1+x_1\right)\exp \left(-{\theta }_1{x}_1\right)\cdots \frac{{\theta }_1^2}{{\theta }_1+1}\left(1+x_k\right)\exp \left(-{\theta }_1{x}_k\right)\\ \cdot \frac{{\theta }_2^2}{{\theta }_2+1}\left(1+x_{k+1}\right)\exp \left(-{\theta }_2{x}_{k+1}\right)\cdots \frac{{\theta }_2^2}{{\theta }_2+1}\left(1+x_n\right)\exp \left(-{\theta }_2{x}_n\right).\\ =\frac{{\theta }_1^{2k}}{{\left({\theta }_1+1\right)}^k}\underset{i=1}{\overset{k}{{\displaystyle \prod }}}\left(1+x_i\right)\exp \left(-{\theta }_1\underset{i=1}{\overset{k}{{\displaystyle \sum }}}{x}_i\right)\frac{{\theta }_2^{2\left(n-k\right)}}{{\left({\theta }_2+1\right)}^{n-k}}\\ \cdot \underset{i=k+1}{\overset{n}{{\displaystyle \prod }}}\left(1+x_i\right)\exp \left(-{\theta }_2\underset{i=k+1}{\overset{n}{{\displaystyle \sum }}}{x}_i\right).\end{array}$

1) 通过Jeffreys提出的用Fisher信息阵来确定 ${\theta }_1,{\theta }_2$ 的无信息先验分布。

样本对数似然函数为：

$\begin{array}{c}l\left({\theta }_1,{\theta }_2|x\right)=2k\ln {\theta }_1-k\ln \left({\theta }_1+1\right)+\underset{i=1}{\overset{k}{{\displaystyle \sum }}}\ln \left(1+x_i\right)-{\theta }_1\underset{i=1}{\overset{k}{{\displaystyle \sum }}}{x}_i+2\left(n-k\right)\ln {\theta }_2\\ -\left(n-k\right)\ln \left({\theta }_2+1\right)+\underset{i=k+1}{\overset{n}{{\displaystyle \sum }}}\ln \left(1+x_i\right)-{\theta }_2\underset{i=k+1}{\overset{n}{{\displaystyle \sum }}}{x}_i.\end{array}$

其中 $x=x_i\left(i=1,2,\cdots ,n\right)$。

通过样本对数似然函数可以求得：

$\frac{\partial l}{\partial {\theta }_1}=\frac{2k}{{\theta }_1}-\frac{k}{{\theta }_1+1}-\underset{i=1}{\overset{k}{{\displaystyle \sum }}}{x}_i,\frac{{\partial }^{2}l}{\partial {\theta }_1^2}=-\frac{2k}{{\theta }_1^2}+\frac{k}{{\left({\theta }_1+1\right)}^2},\frac{{\partial }^{2}l}{\partial {\theta }_1\partial {\theta }_2}=0.$

$\frac{\partial l}{\partial {\theta }_2}=\frac{2\left(n-k\right)}{{\theta }_2}-\frac{n-k}{{\theta }_2+1}-\underset{i=k+1}{\overset{n}{{\displaystyle \sum }}}{x}_i,\frac{{\partial }^{2}l}{\partial {\theta }_2^2}=-\frac{2\left(n-k\right)}{{\theta }_2^2}+\frac{n-k}{{\left({\theta }_2+1\right)}^2},\frac{{\partial }^{2}l}{\partial {\theta }_2\partial {\theta }_1}=0.$

进而得到 ${\theta }_1,{\theta }_2$ 的无信息先验矩阵为：

$I\left({\theta }_1,{\theta }_2\right)=E^{x|{\theta }_1,{\theta }_2}\left(-\frac{{\partial }^{2}l}{\partial {\theta }_i\partial {\theta }_j}\right)=\left(\begin{array}{cc}\frac{2k}{{\theta }_1^2}-\frac{k}{{\left({\theta }_1+1\right)}^2}& 0\\ 0& \frac{2\left(n-k\right)}{{\theta }_2^2}-\frac{n-k}{{\left({\theta }_2+1\right)}^2}\end{array}\right).$

其中 $i,j=1,2$。

故 ${\theta }_1,{\theta }_2$ 的无信息先验分布为：

$\begin{array}{c}\pi \left({\theta }_1,{\theta }_2\right)={\left[\det I\left({\theta }_1,{\theta }_2\right)\right]}^{\frac{1}{2}}={\left[\left(\frac{2k}{{\theta }_1^2}-\frac{k}{{\left({\theta }_1+1\right)}^2}\right)\left(\frac{2\left(n-k\right)}{{\theta }_2^2}-\frac{n-k}{{\left({\theta }_2+1\right)}^2}\right)\right]}^{\frac{1}{2}}\\ ={\left[\frac{k\left(n-k\right)\left[2{\left({\theta }_1+1\right)}^2-{\theta }_1^2\right]\left[2{\left({\theta }_2+1\right)}^2-{\theta }_2^2\right]}{{\theta }_1^2{\left({\theta }_1+1\right)}^2{\theta }_2^2{\left({\theta }_2+1\right)}^2}\right]}^{\frac{1}{2}}.\end{array}$

由贝叶斯公式求得 $k,{\theta }_1,{\theta }_2$ 的联合后验分布为：

$\begin{array}{l}\pi \left(k,{\theta }_1,{\theta }_2|x\right)\propto L\left(k,{\theta }_1,{\theta }_2|x\right)\pi \left(k\right)\pi \left({\theta }_1,{\theta }_2\right)\\ =\frac{{\theta }_1^{2k}}{{\left({\theta }_1+1\right)}^k}\underset{i=1}{\overset{k}{{\displaystyle \prod }}}\left(1+x_i\right)\exp \left(-{\theta }_1\underset{i=1}{\overset{k}{{\displaystyle \sum }}}{x}_i\right)\cdot \frac{{\theta }_2^{2\left(n-k\right)}}{{\left({\theta }_2+1\right)}^{n-k}}\cdot \underset{i=k+1}{\overset{n}{{\displaystyle \prod }}}\left(1+x_i\right)\exp \left(-{\theta }_2\underset{i=k+1}{\overset{n}{{\displaystyle \sum }}}{x}_i\right)\\ \cdot \frac{1}{n-1}\cdot {\left[\frac{k\left(n-k\right)\left[2{\left({\theta }_1+1\right)}^2-{\theta }_1^2\right]\left[2{\left({\theta }_2+1\right)}^2-{\theta }_2^2\right]}{{\theta }_1^2{\left({\theta }_1+1\right)}^2{\theta }_2^2{\left({\theta }_2+1\right)}^2}\right]}^{\frac{1}{2}}.\end{array}$

各参数满条件分布为：

$f\left({\theta }_1|{\theta }_2,k,x_i\right)\propto \frac{{\theta }_1^{2k}}{{\left({\theta }_1+1\right)}^k}\exp \left(-{\theta }_1\underset{i=1}{\overset{k}{{\displaystyle \sum }}}{x}_i\right)\cdot {\left[\frac{\left[2{\left({\theta }_1+1\right)}^2-{\theta }_1^2\right]}{{\theta }_1^2{\left({\theta }_1+1\right)}^2}\right]}^{\frac{1}{2}};$

$f\left({\theta }_2|{\theta }_1,k,x_i\right)\propto \frac{{\theta }_2^{2\left(n-k\right)}}{{\left({\theta }_2+1\right)}^{n-k}}\exp \left(-{\theta }_2\underset{i=k+1}{\overset{n}{{\displaystyle \sum }}}{x}_i\right)\cdot {\left[\frac{\left[2{\left({\theta }_2+1\right)}^2-{\theta }_2^2\right]}{{\theta }_2^2{\left({\theta }_2+1\right)}^2}\right]}^{\frac{1}{2}};$

$f\left(k|{\theta }_1,{\theta }_2,x_i\right)\propto \frac{{\theta }_1^{2k}}{{\left({\theta }_1+1\right)}^k}\underset{i=1}{\overset{k}{{\displaystyle \prod }}}\left(1+x_i\right)\exp \left(-{\theta }_1\underset{i=1}{\overset{k}{{\displaystyle \sum }}}{x}_i\right)\cdot \frac{{\theta }_2^{2\left(n-k\right)}}{{\left({\theta }_2+1\right)}^{n-k}}$

$\underset{i=k+1}{\overset{n}{{\displaystyle \prod }}}\left(1+x_i\right)\exp \left(-{\theta }_2\underset{i=k+1}{\overset{n}{{\displaystyle \sum }}}{x}_i\right)\cdot {\left[k\left(n-k\right)\right]}^{\frac{1}{2}}.$

比较选取均匀分布作为先验分布来说，Jeffreys提出的用Fisher信息阵来确定 ${\theta }_1,{\theta }_2$ 的无信息先验分布在单调变换中具有不变性，能够保证不论采取什么样的参数化方法，它们的先验分布始终是互通的，从而后验分布也是互通的。

2) ${\theta }_1,{\theta }_2$ 的先验分布为伽玛分布 $\left(Ga\left(b_i,c_i\right),i=1,2\right)$。

$\{\begin{array}{l}\pi \left({\theta }_1\right)=Ga\left({\theta }_1|b_1,c_1\right)=\frac{c_1^{b_1}}{\Gamma \left(b_1\right)}{\theta }_1^{b_1-1}{\text{e}}^{-c_1{\theta }_1},{\theta }_1>0,b_1>0,c_1>0.\\ \pi \left({\theta }_2\right)=Ga\left({\theta }_2|b_2,c_2\right)=\frac{c_2^{b_2}}{\Gamma \left(b_2\right)}{\theta }_2^{b_2-1}{\text{e}}^{-c_2{\theta }_2},{\theta }_2>0,b_2>0,c_2>0.\end{array}$

且 $k,{\theta }_1,{\theta }_2$ 相互独立，由贝叶斯公式得 $k,{\theta }_1,{\theta }_2$ 的联合后验密度为：

$\begin{array}{c}\pi \left(k,{\theta }_1,{\theta }_2|x\right)\propto L\left(k,{\theta }_1,{\theta }_2|x\right)\pi \left(k\right)\pi \left({\theta }_1\right)\pi \left({\theta }_2\right)\\ =\frac{1}{n-1}\cdot \frac{{\theta }_1^{2k}}{{\left({\theta }_1+1\right)}^k}\underset{i=1}{\overset{k}{{\displaystyle \prod }}}\left(1+x_i\right)\cdot \exp \left(-{\theta }_1\underset{i=1}{\overset{k}{{\displaystyle \sum }}}{x}_i\right)\frac{{\theta }_2^{2\left(n-k\right)}}{{\left({\theta }_2+1\right)}^{n-k}}\underset{i=k+1}{\overset{n}{{\displaystyle \prod }}}\left(1+x_i\right)\\ \cdot \exp \left(-{\theta }_2\underset{i=k+1}{\overset{n}{{\displaystyle \sum }}}{x}_i\right)\frac{c_1^{b_1}}{\Gamma \left(b_1\right)}{\theta }_1^{b_1-1}{\text{e}}^{-c_1{\theta }_1}\cdot \frac{c_2^{b_2}}{\Gamma \left(b_2\right)}{\theta }_2^{b_2-1}{\text{e}}^{-c_2{\theta }_2}.\end{array}$

各参数满条件分布为：

$f\left({\theta }_1|{\theta }_2,k,x_i\right)\propto \frac{{\theta }_1^{2k}}{{\left({\theta }_1+1\right)}^k}\exp \left(-{\theta }_1\underset{i=1}{\overset{k}{{\displaystyle \sum }}}{x}_i\right)\cdot {\theta }_1^{b_1-1}{\text{e}}^{-c_1{\theta }_1};$

$f\left({\theta }_2|{\theta }_1,k,x_i\right)\propto \frac{{\theta }_2^{2\left(n-k\right)}}{{\left({\theta }_2+1\right)}^{n-k}}\exp \left(-{\theta }_2\underset{i=k+1}{\overset{n}{{\displaystyle \sum }}}{x}_i\right)\cdot {\theta }_2^{b_2-1}{\text{e}}^{-c_2{\theta }_2};$

$f\left(k|{\theta }_1,{\theta }_2,x_i\right)\propto \frac{{\theta }_1^{2k}}{{\left({\theta }_1+1\right)}^k}\underset{i=1}{\overset{k}{{\displaystyle \prod }}}\left(1+x_i\right)\exp \left(-{\theta }_1\underset{i=1}{\overset{k}{{\displaystyle \sum }}}{x}_i\right)\frac{{\theta }_2^{2\left(n-k\right)}}{{\left({\theta }_2+1\right)}^{n-k}}\underset{i=k+1}{\overset{n}{{\displaystyle \prod }}}\left(1+x_i\right)\exp \left(-{\theta }_2\underset{i=k+1}{\overset{n}{{\displaystyle \sum }}}{x}_i\right).$

4. 随机模拟

在随机模拟过程中，考虑到参数 ${\theta }_1,{\theta }_2,k$ 的满条件分布比较复杂，因此选用M-H算法对各参数的满条件分布进行抽样。接下来介绍Markov Chain Monte Carlo (MCMC)算法的几个具体步骤：

设初始点 ${\alpha }^{\left(0\right)}=\left(k^{\left(0\right)},{\theta }_1^{\left(0\right)},{\theta }_2^{\left(0\right)}\right)$ 经过迭代后第 $t-1$ 次迭代值为 ${\alpha }^{\left(t-1\right)}=\left(k^{\left(t-1\right)},{\theta }_1^{\left(t-1\right)},{\theta }_2^{\left(t-1\right)}\right)$，则第t次迭代步骤如下：

1) ${\theta }_1^{\left(t\right)}~f\left({\theta }_1|{\theta }_2,k,x_i\right)$，选取建议分布 $q\left({\theta }_1^{\left(t-1\right)},{{\theta }^{\prime }}_1\right)$ 为均匀分布，并从中随机抽取 ${{\theta }^{\prime }}_1$，令 $r\left({\theta }_1^{\left(t-1\right)},{{\theta }^{\prime }}_1\right)=\min \left\{\frac{\pi \left({{\theta }^{\prime }}_1|\cdot \right)}{\pi \left({\theta }_1^{\left(t-1\right)}|\cdot \right)},1\right\}$，若随机数 $u\le r\left({\theta }_1^{\left(t-1\right)},{{\theta }^{\prime }}_1\right)$，则 ${\theta }_1^{\left(t\right)}={{\theta }^{\prime }}_1$，否则 ${\theta }_1^{\left(t\right)}={\theta }_1^{\left(t-1\right)}$ ；

2) ${\theta }_2^{\left(t\right)}~f\left({\theta }_2|{\theta }_1,k,x_i\right)$，获取 ${\theta }_2^{\left(t\right)}$ 与1)类似；

3) $k^{\left(t\right)}~f\left(k|{\theta }_1,{\theta }_2,x_i\right)$，选取建议分布 $q\left(k^{\left(t-1\right)},k^{\prime }\right)$ 为取值 $0,1,\cdots ,n-1$ 的离散型均匀分布，并从中随机抽取 $k^{\prime }$，令 $r\left(k^{\left(t-1\right)},k^{\prime }\right)=\min \left\{\frac{\pi \left(k^{\prime }|\cdot \right)}{\pi \left(k^{\left(t-1\right)}|\cdot \right)},1\right\}$，若随机数 $u\le r\left(k^{\left(t-1\right)},k^{\prime }\right)$，则 $k^{\left(t\right)}=k^{\prime }$，否则 $k^{\left(t\right)}=k^{\left(t-1\right)}$。

设 $\left(k^{\left(j\right)},{\theta }_1^{\left(j\right)},{\theta }_2^{\left(j\right)}\right),j=1,2,\cdots ,B,\cdots ,M$ 为迭代M次所得的Gibbs样本，若B次后迭代逐渐收敛，则将后 $M-B$ $M-B$ 个迭代的均值作为参数 $k,{\theta }_1,{\theta }_2$ 的估计值，

$\hat{k}=\frac{1}{M-B}\underset{t=B+1}{\overset{M}{{\displaystyle \sum }}}{k}^{\left(t\right)},{\hat{\theta }}_i=\frac{1}{M-B}\underset{t=B+1}{\overset{M}{{\displaystyle \sum }}}{\theta }_i{}^{\left(t\right)},i=1,2.$

取 $n=200$ 个样本，参数 $\left(k,{\theta }_1,{\theta }_2\right)$ 的真实值取 $\left(100,3,8\right)$，此时Lindley分布的模型为：

$X_i~\{\begin{array}{l}\frac{9}{4}\left(1+x\right)\exp \left(-3x\right),i=1,2,\cdots ,100\\ \frac{64}{9}\left(1+x\right)\exp \left(-8x\right),i=101,102,\cdots ,200\end{array}$

利用各参数的满条件分布，运用R软件进行MCMC模拟。为确保参数的收敛性，先进行10,000次的预迭代，再进行20,000次迭代。结果如下所示：

1) 当 ${\theta }_1,{\theta }_2$ 选取无信息先验分布时：

Table 1. Bayesian estimation of parameters k , θ 1 , θ 2 under uninformative prior distribution

表1. 无信息先验分布下参数 $k,{\theta }_1,{\theta }_2$ 的贝叶斯估计

Figure 1. The iteration trajectory of parameter k

图1. 参数k的迭代轨迹

Figure 2. The iteration trajectory of parameter k

图2. 参数k的迭代轨迹

2) 当 ${\theta }_1,{\theta }_2$ 选取伽玛先验分布时： ${\theta }_1~Ga\left(7,3\right)$， ${\theta }_2~Ga\left(9,2\right)$。

Table 2. Bayesian estimation of parameters k , θ 1 , θ 2 under conjugate prior distribution

表2. 共轭先验分布下参数 $k,{\theta }_1,{\theta }_2$ 的贝叶斯估计

Figure 3. Two iteration trajectories of the parameter k

图3. 参数k的两条迭代图

Figure 4. Two iteration trajectories of the parameter k

图4. 参数k的两条迭代轨迹

结果分析：由表1和表2知，当参数选取不同的先验分布后再进行随机模拟，得到各参数估计值与真实值的MC误差均不超过2%，因此各参数的估计值在较高水平上是有效的；各参数置信水平0.95的置信区间[2.5%分位数，97.5%分位数]较窄，说明区间估计效果良好；图1，图2是变点k的抽样迭代轨迹，可以根据图上信息判断样本是否收敛。两张图上显示出抽样基本都在变点附近波动，具有一定的规律性；此外，由图3，图4看出k的两条Markov链趋于重合，具有较好的收敛性。综上可得，Lindley分布的参数和变点估计可由MCMC算法得到较为理想的效果，可用该方法解决Lindley分布的变点问题。

基金项目

国家自然科学基金项目(11801488)；新疆师范大学教学研究与改革项目(SDJG2020-30)；新疆师范大学科研发展专项项目(XJNUZX202001)。

NOTES

^*第一作者。

^#通讯作者。