冲击模型是可靠性数学理论中的重要研究内容之一, 经典的冲击模型以冲击次数为失效机制, 而δ冲击模型是以冲击间隔为失效机制的.自从Li提出δ冲击模型之后^[1], 有关δ方面的文章日益增多.

δ冲击模型是指相邻两次冲击的时间间隔小于事先给定的参数(或称为阈值)δ时, 系统发生故障. Li讨论了δ冲击模型的概率分布与寿命性质^[2]; 李泽慧等针对一般冲击模型, 研究Bayes方法处理无失效数据的问题, 分别选择均匀分布和Beta分布作为先验分布, 用Bayes方法和多层Bayes方法得到了参数$\delta_1$和$\delta_2$的估计^[3]; 马明讨论了δ冲击模型寿命分布中多重积分的计算问题, 并研究了由此积分引出的M函数的性质^[4]; 冶建华等研究了时间点服从$0 \sim 1$分布的离散型δ冲击模型系统寿命的概率分布和期望性质^[5].有关δ冲击模型的进一步发展, 可参见文献[6-10].

最近Ma扩展了δ冲击模型, 提出了截断δ冲击模型^[11].截断δ冲击模型是指相邻两次冲击的时间间隔大于事先给定的参数δ时, 系统失效.张攀等研究了时间点服从0--1分布的离散型截断δ冲击模型中系统寿命的分布和期望^[12].郑莹等讨论了自激滤过的泊松过程的二阶矩, 并将所得结论应用于截断δ冲击模型的标值过程和客户寿命价值, 得到了截断δ冲击模型的标值过程的二阶矩和客户寿命价值的二阶矩^[13]. Serkan Eryilmaz研究了冲击时间间隔服从更新过程的截断δ冲击模型, 得到寿命的前1、2阶矩, 进一步得到了更新间隔为均匀分布情形时的生存函数^[14].

前人文章中对截断δ冲击模型的研究主要关注于寿命分布的计算, 而对参数估计问题涉足比较少.本文在三种样本轨道数据的情形下讨论了截断δ冲击模型的完全寿命参数估计问题.

假设系统遭受的外部冲击时间间隔服从参数为$\lambda>0$的指数分布, 而且任意两次的冲击时间间隔是相互独立的.设$X_n (n=1, 2, \cdots)$为第$(n-1)$次和第n次冲击的时间间隔, $t_n$表示第n次冲击的到达时间, $\delta>0$是一个实数, 则截断δ冲击模型的寿命可以描述为:

$\begin{eqnarray} \overline{T}=\sum\limits_{n=0}^{\overline{N}}X_n +\delta =t_{\overline{N}} +\delta, \nonumber \end{eqnarray}$

其中$\overline{N}=\inf\{{n: {t_{n+1}-t_n} \geq \delta, n \geq 0}\}.$

对于寿命$\overline{T}$有下面的引理:

引理 ^[11] 截断$\delta$冲击模型寿命的期望和二阶矩分别为

$E[\overline{T}]={\frac{1}{\lambda}}(e^{{\lambda}{\delta}}-1)\nonumber \:\mbox{和}\: E[\overline{T}^2]={\frac{2}{\lambda}}{e^{\lambda\delta}}{[\frac{1}{\lambda}(e^{\lambda\delta}-1)-\delta]}.\nonumber$

在截断δ冲击模型中, $\delta$和$\lambda$是两个重要的参数.为了估计$\delta$和$\lambda$, 我们假设观察出的样本轨道数据有下面三种形式, 并且是完全寿命数据, 即全部观测到寿命终止为止.

假设$A_{1}$:观察到一列冲击到达时间数据$t_{0}, t_{1}, t_{2}, \cdots, t_{n}, n\geq 1$且$t_{n}$是寿命终止前最后一次冲击时刻(其中$t_{0}=0$).

假设$A_{2}$:观测到寿命终止时一共冲击次数是n, $n\geq1$.

假设$A_{3}$:观察到多组寿命结束时间$T_{1}, T_{2}, \cdots, T_{m}, m \geq 1$.

注   1)  在假设$A_{1}$下, 意味着$t_{n+1}-t_{n}>\delta$, 且若$n\geq1, t_{i+1}-t_{i}<\delta, i=0, \cdots, n-1$, 但$t_{n+1}$不能观测到.

2)  此外, 若允许n=0, 即只有一个观测值$t_{0}=0$, 此时是不可识别的, 因为我们不清楚是纯粹没有冲击发生还是有冲击到达而被δ截断了.

定理1 在假设$A_{1}$下, 若$\delta$已知, 则$\lambda$的极大似然估计为

$\begin{equation} \label{EL} \widehat{\lambda}=\frac{n}{t_{n}+\delta}. \end{equation}$

(2.1)

证 在假设$A_{1}$下, 由于时间间隔服从指数分布, 并且在$t_{n}$后附加了一段大于$\delta$的时间, 所以样本的似然函数为

$L(\lambda ) = \left[ {\prod\limits_{i = 1}^n {(\lambda {e^{ - \lambda ({t_i} - {t_{i - 1}})}})} } \right]{e^{ - \lambda \delta }} = {\lambda ^n}{e^{ - \lambda {t_n}}}{e^{ - \lambda \delta }}.$

对似然函数取对数得

$\begin{equation} \label{LNLL} \ln(L(\lambda))=n\ln\lambda-{\lambda}t_n-\lambda\delta. \end{equation}$

(2.2)

对式$(2.2)$求关于$\lambda$的导数, 并令其等于0得$\widehat{\lambda}=\frac{n}{t_{n}+\delta}.$

定理得证.

定理2  在假设$A_{1}$下, 若$\delta$与$\lambda$均未知, 则$\delta$与$\lambda$的极大似然估计分别是

$\widehat{\delta}=\max\{t_i-t_{i-1}; i=1, \cdots, n\}\:\mbox{和}\: \widehat{\lambda}=\frac{n}{t_n+\max\{t_i-t_{i-1}; i=1, \cdots, n\}}\:, \:n>0.$

其中, $t_0=0$.

证  对式$(2.2)$关于$\delta$与$\lambda$分别取导数, 得

$\begin{eqnarray} \label{PLL} \frac{\partial\ln{L}}{\partial\lambda}=\frac{n}{\lambda}-t_n-\delta, \end{eqnarray}$

(2.3)

$\begin{eqnarray} \label{PLD} \frac{\partial\ln{L}}{\partial\delta}={-\lambda}. \end{eqnarray}$

(2.4)

由于式$(2.4)$小于0, 说明$\delta$取最小值时似然函数最大, 但注意到样本数据中的每个时间间隔都小于$\delta$, 所以得$\widehat{\delta}=\max\{t_i-t_{i-1}; i=1, \cdots, n\}.$令式$(2.3)$等于0, 将$\widehat{\delta}$代入得

$\widehat{\lambda}=\frac{n}{t_n+\max\{t_i-t_{i-1}; i=1, \cdots, n\}}.$

注 实际应用中, 取$\delta+\varepsilon$, 其中$\varepsilon>0$是微小的修正量.

定理3 设$\widehat \delta $是假设$A_{1}$下$\delta$的极大似然估计, 则

$E[\widehat{\delta}]=\delta-\frac{1}{\lambda}(1-e^{-\lambda\delta}).$

证 对$\overline{N}$的任一实现n, 记$X_{i}=t_{i}-t_{i-1}, i=1, {\cdots}, n+1$ (其中$X_{0}=0$), 有

$ E[\widehat{\delta}|\overline{N}=n]=E[\max\{X_{i}, i=1, {\cdots}, n\}|X_{i}＜\delta; i=1, {\cdots}, n, X_{n+1}＞\delta].$

注意到$\max\{X_{i}, i=1, {\cdots}, n\}＜t\Longleftrightarrow\{X_{i}＜t, i=1, {\cdots}, n\}, $并且当$t\geq\delta$时,

$P(X_{i}＜t, i=1, {\cdots}, n|X_{i}＜\delta; i=1, {\cdots}, n, X_{n+1}>\delta)=1, $

当$0<t<\delta$时,

$ P(X_{i}＜t, i=1, {\cdots}, n|X_{i}＜\delta; i=1, {\cdots}, n, X_{n+1}＞\delta) \\ =\frac{P(X_{i}＜t, i=1, {\cdots}, n, X_{i}＜\delta; i=1, {\cdots}, n, X_{n+1}＞\delta)}{P(X_{i}＜\delta; i=1, {\cdots}, n, X_{n+1}＞\delta)}\\ =\frac{P(X_{i}＜t, i=1, {\cdots}, n, X_{n+1}＞\delta)}{P(X_{i}＜\delta; i=1, {\cdots}, n, X_{n+1}＞\delta)}=\frac{\prod\limits_{i=1}^{n}(1-e^{-\lambda t})}{\prod\limits_{i=1}^{n}(1-e^{-\lambda\delta})}=\frac{(1-e^{-\lambda t})^n}{(1-e^{-\lambda\delta})^n}, $

则

$ E[\widehat{\delta}|\overline{N}=n]=\int_0^{+\infty}P(\widehat{\delta}>t|\overline{N}=n)\text{d}t=\int_0^{\delta}[1-\frac{(1-e^{-\lambda t})^n}{(1-e^{-\lambda\delta})^n}]\text{d}t=\delta-\int_0^{\delta}{\frac{(1-e^{-\lambda t})^n}{(1-e^{-\lambda\delta})^n}}\text{d}t.$

由全期望公式, 得

$E[\widehat{\delta}]=\mathop \sum \limits_{n = 0}^\infty [\delta-\int_0^{\delta}{\frac{(1-e^{{-\lambda}t})^n}{(1-e^{-\lambda\delta})^n}}\text{d}t](1-e^{-\lambda\delta})^ne^{-\lambda\delta}\\ ={\delta}e^{-\lambda\delta}\mathop \sum \limits_{n = 0}^\infty (1-e^{-\lambda\delta})^n-e^{-\lambda\delta}\mathop \sum \limits_{n = 0}^\infty \int_0^{\delta}(1-e^{-\lambda{t}})^n\text{d}t\\ =\delta-e^{-\lambda\delta}\int_0^{\delta}{\mathop \sum \limits_{n = 0}^\infty }(1-e^{-\lambda{t}})^n\text{d}t=\delta-\frac{1}{\lambda}(1-e^{-\lambda\delta}).$

定理得证.

定理3表明, 在假设$A_{1}$下, $\delta$的极大似然估计是有偏估计.我们可以由定理3结果构造$\delta$的一个无偏估计.

定理4  在假设$A_{1}$下, 若$\lambda$已知, 则

$\max\{t_{i}-t_{i-1}; i=1, {\cdots}, \overline{N}\}+\frac{1}{\lambda}I_{\{\overline{N}>0\}}$

是$\delta$的一个无偏估计.

证 由定理3有

$E[\max\{t_{i}-t_{i-1};i=1, {\cdots}, \overline{N}\}]=\delta-\frac{1}{\lambda}(1-e^{-\lambda\delta}).$

所以我们只需证明$E[\frac{1}{\lambda}I_{\{\overline{N}>0\}}]=\frac{1}{\lambda}(1-e^{-\lambda\delta}), $而

$ E[I_{\{\overline{N}>0\}}]=P(\overline{N}>0)=1-P(\overline{N}=0)=1-P(X_{i}>\delta)=1-e^{(-\lambda\delta)}, $

定理得证.

在假设$A_{2}$下, 我们有下面的估计.

定理5  在假设$A_{2}$下,

$1)$若$\delta$已知, 则$\lambda$的极大似然估计是$\widehat{\lambda}=\frac{\ln(n+1)}{\delta}$;

$2)$若$\lambda$已知, 则$\delta$的极大似然估计是$\widehat{\delta}=\frac{\ln(n+1)}{\lambda}$;

$3)$若$\delta$与$\lambda$都未知, 则$\lambda\delta$的极大似然估计是$\widehat{\lambda\delta}=\ln(n+1)$.

证 我们先求似然函数, 寿命结束时一共冲击n次, 则

$L(\delta, \lambda)=P(\overline{N}=n)=P\{\bigcap\limits_{i=1}^n(X_{i}＜\delta)(X_{n+1}＞\delta)\}.$

所以

$L(\delta, \lambda)=(1-e^{-\lambda\delta})^ne^{-\lambda\delta}.$

分别取对数求偏导并令等于0得

$\begin{equation} \label{eq:1} \left\{ \begin{aligned} \frac{{\partial}L}{\partial\delta}=n\frac{{\lambda}e^{-\lambda\delta}}{1-e^{-\lambda\delta}}-\lambda=0, \nonumber\\ \frac{{\partial}L}{\partial\lambda}=n\frac{{\delta}e^{-\lambda\delta}}{1-e^{-\lambda\delta}}-\delta=0, \nonumber\\ \frac{{\partial}L}{\partial\lambda\delta}=n\frac{e^{-\lambda\delta}}{1-e^{-\lambda\delta}}-1=0.\nonumber \end{aligned} \right. \end{equation}$

分别解方程得证.

注  若将$\widehat{\delta}=\frac{\ln(n+1)}{\lambda}$直接代入$E[\overline{T}]=\frac{1}{\lambda}(e^{\lambda\delta}-1)$可得假设$A_{2}$下, 平均寿命的一个估计式$\widehat{E[\overline{T}]}=\frac{n}{\lambda}$, 这恰如参数为n和$\lambda$的伽玛分布的期望, 这个估计有一个直观的解释, 即平均寿命等于平均冲击间隔乘以总冲击次数.此外, 由于$\overline{N}$服从几何分布, $E[\overline{N}]=e^{\lambda\delta}-1$, 所以$E[\widehat{E[\overline{T}]}]=E[\frac{\overline{N}}{\lambda}]=\frac{1}{\lambda}(e^{\lambda\delta}-1)$, 因此$\frac{\overline{N}}{\lambda}$也是$E[\overline{T}]$的一个无偏估计.

在假设$A_{3}$下, 我们有下面的估计.

定理6  在假设$A_{3}$下, 记$\overline{X}=\frac{1}{m}\sum\limits_{i=1}^{m}T_{i}, \overline{X^2}=\frac{1}{m}\sum\limits_{i=1}^{m}{T_{i}}^2$, 则

$1)$若$\lambda$已知, 则$\delta$的矩估计量为$\widehat{\delta}=\frac{\ln(\lambda\overline{X}+1)}{\lambda}$;

$2)$若$\delta$已知, 则$\lambda$的矩估计量为$\widehat{\lambda}=\frac{2(\overline{X}-\delta)}{\overline{X^2}-2\overline{X}(\overline{X}-\delta)}$.

证 令

$\left\{ \begin{array}{l} \bar X = E[\bar T],\\ \overline {{X^2}} = E[{{\bar T}^2}]. \end{array} \right.$

由引理结论得证.

下面我们应用截断δ冲击模型对客户关系营销系统中的客户到达率进行估计.

考虑如下泊松营销系统

${\rm{US /CM /C /M /i}}{\rm{.id / \mathsf{ δ} }},$

其中/US/意味着交易类型是不依订单可任意时刻到达的; /CM/表明客户行为是迁移型的; /C/表明每期的营销支出是一个常数, 不妨设为M; /M/说明响应是一个泊松过程, 满足无记忆性, 设其到达率为$\lambda$; /i. i. d/表示每次净贡献是一个独立同分布的序列; $/\delta/$意味着中止关系采用的策略是δ策略. Ma et al. (2008) 建立了营销系统的数学框架, 并讨论了马尔科夫营销系统的最优营销决策^[15], 郑莹等(2013) 讨论了泊松营销系统的二阶矩^[13], 有关关系营销系统的进一步研究可参见文献[13]和[15].

假设再营销仅影响顾客的到达率, 即$\lambda$是M的一个函数$\lambda(M)$, 称为到达率, 常用的一个到达率函数是

$\begin{eqnarray} \label{LFM} \lambda(M)={\rm arr}(1-e^{-{\theta}M}), \end{eqnarray}$

(3.1)

式$(3.1)$中, $\theta$表示公司目前保留的水平值, arr表示客户到达最大可能人数. Blattberg et al.(1996) 使用式$(3.1)$计算了客户资产(Customer Equity)\ucite{RAJ}, Ma et al.(2008) 使用式$(3.1)$计算了$\delta$和M的最优值^[15].

在泊松关系营销系统中, 参数估计主要的任务是在不同的营销支出M下, 根据客户的到达规律估计其到达率$\lambda$, 估计的目的是为了进行营销决策.估计中需要的样本数据是不同的营销支出下, 客户的响应行为样本轨道.但营销有其特殊性, 即客户的响应行为是不可再现的, 这是由于, 若得到了一个样本数据, 说明该客户关系的寿命已结束了.所以, 在收集样本数据时, 一般是选择与当前客户各种社会自然特征相近的典型客户.

假设我们得到的数据是不同的营销支出${M}$下, 相应典型客户的所有相应时间(完全寿命), 如表 1所示:

根据文献[11], 客户响应行为就是一个典型的截断δ冲击模型, 其中客户响应时间相当于冲击到达时刻, 表 1中每一行响应时间就相当于假设$A_{1}$的形式, 所以由$(2.1)$式很容易估计出客户的到达率, 即

$\begin{eqnarray} \label{LIE} \widehat{\lambda_{i}}=\frac{n_{i}}{t_{i, n_{i}}+\delta}, \; n=1, \cdots, m. \end{eqnarray}$

(3.2)

为了检验估计的效果, 我们按照如下路线进行估计:在不同的营销支出下, 选择适当的参数$\theta$和arr, 根据式$(3.1)$计算出相应的到达率$\lambda$, 应用这些$\lambda$生成客户的响应时间轨道(按照泊松营销系统的定义, 响应时间间隔服从指数分布), 从而得到表 1的具体数值, 然后根据这些数值代入式$(3.2)$得到估计值$\widehat{\lambda}$.

我们选取$\theta$= 0.03, arr = 0.75, M按照步长为0.5在$1\sim10$之间取值, 运用MATLAB编程实现, 为了减少随机误差, 程序中使用了批均值法, 结果见表 2.

表 2中的第2列是由式$(3.1)$得到的$\lambda$真值, 第3列至第7列是基于不同的$\delta$下所截断的响应轨道得到的$\lambda$的估计值. 图 1表示了$\lambda$的估计值与$\lambda$真值之间的误差.从图 1可看出, 估计值围绕真值曲线进行波动, 误差基本没有超过±0.1.此外, 随着$\delta$的增加, 误差反而减小.

表 2中, 每一行的估计值是针对同一个$\lambda$进行估计的, 为了更进一步直观考察估计效果, 对每一行估计值计算平均值$\bar{\lambda}$, 将$\lambda$和$\bar{\lambda}$关于${M}$变化曲线分别绘于图 2中, 从图 2可看出平均后的最大误差不超过0.04 (在${M}$=3处).