模型推导:BG/NBD(预测用户生命周期(CLV)模型)

CLV(Customer Lifetime Value)指的是客户生命周期价值,用以衡量客户在一段时间内对企业有多大的价值。企业对每个用户的流失与否、在未来时间是否会再次购买,还会再购买多少次才会流失等问题感兴趣,本文中的BG/NBD模型就是用来解决这样一系列问题的。

本文的模型数学推理均参考自:
https://www.brucehardie.com/notes/039/bgnbd_derivation__2019-11-06.pdf

模型假设

BG/NBD模型包含以下的五个假设:
1、用户在活跃状态下,一个用户在时间段t内完成的交易数量服从交易率为λ\lambdaλ的泊松过程,这相当于假设了交易与交易之间的时间随交易率λ\lambdaλ呈指数分布

f(tj∣tj−1;λ)=λe−λ(tj−tj−1)f(t_j|t_{j-1};\lambda) = \lambda e^{-\lambda(t_j-t_{j-1})}f(tjtj1;λ)=λeλ(tjtj1) , t_j>t_{j-1}>=0

2、用户的交易率λ\lambdaλ服从如下的gamma分布

f(λ∣r,α)=αrλr−1e−λαΓ(r),λ>0f(\lambda | r,\alpha)=\frac{\alpha^r\lambda ^{r-1}e^{-\lambda\alpha}}{\Gamma(r)} , \lambda>0f(λr,α)=Γ(r)αrλr1eλα,λ>0

其中Γ(r)=∫0+∞tr−1e−tdt\Gamma(r)=\int_0^{+\infty}t^{r-1}e^{-t}dtΓ(r)=0+tr1etdt,是gamma函数

3、每个用户在第j次交易完之后流失的概率服从参数为p的二项分布,p为发生任何交易后用户流失的概率

P(在第j次交易后不再活跃)=p(1−p)j−1,j=1,2,3,...P(在第j次交易后不再活跃)=p(1-p)^{j-1},j=1,2,3,...P(在第j次交易后不再活跃)=p(1p)j1,j=1,2,3,...

4、每个用户的流失率p服从形状参数为a,b的beta分布

f(p∣a,b)=pa−1(1−p)b−1B(a,b),0<=p<=1f(p|a,b)=\frac{p^{a-1}(1-p)^{b-1}}{B(a,b)},0<=p<=1f(pa,b)=B(a,b)pa1(1p)b10<=p<=1

其中,B(a,b)是贝塔函数,B(a,b)=Γ(a)Γ(b)Γ(a+b)=∫01xa−1(1−x)b−1dxB(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}=\int_0^1x^{a-1}\left(1-x\right)^{b-1}dxB(a,b)=Γ(a+b)Γ(a)Γ(b)=01xa1(1x)b1dx

5、交易率λ\lambdaλ和流失率p互相独立

对于某个特定顾客的推导(已知λ\lambdaλ和p)

1、似然函数的推导

t1t_1t1时刻发生第一次交易的概率服从标准指数分布,等于λe−λt1\lambda e^{-\lambda t_1}λeλt1
t2t_2t2时刻发生第二次交易的概率为t1t_1t1时刻依然活跃的概率乘以标准指数分布,等于λ(1−p)λe−λt2−t1\lambda (1-p)\lambda e^{-\lambda t_2-t_1}λ(1p)λeλt2t1
以此类推,
txt_xtx时刻发生第x次交易的概率为txt_xtx时刻依然活跃的概率乘以标准指数分布,等于λ(1−p)λe−λtx−tx−1\lambda (1-p)\lambda e^{-\lambda t_x-t_{x-1}}λ(1p)λeλtxtx1
而在(tx,T](t_x,T](tx,T]时间段内,没有任何交易的概率,则是客户已经不再活跃额度概率加上他活跃但是不进行任何购买的概率,也就是p+(1−p)e−λT−txp+(1-p)e^{-\lambda T-t_x}p+(1p)eλTtx
注意后面那一项是因为泊松过程发生k次事件的概率

P(N(T)−N(tx)=k)=e−λ(T−tx)(λ(T−tx))kk!P(N(T)-N(t_x)=k)=\frac{e^{-\lambda (T-t_x)}(\lambda (T-t_x))^k}{k!}P(N(T)N(tx)=k)=k!eλ(Ttx)(λ(Ttx))k 令k等于0得到,再乘以在t_x时刻依然活跃的概率(1-p)得到的。
因此,似然函数可以写作

L(λ,p∣t1,t2,...,tx,T)=(1−p)x−1λxe−λtx∗(p+(1−p)e−λT−tx)=p(1−p)x−1λxe−λtx+(1−p)xλxe−λL(\lambda ,p|t_1,t_2,...,t_x,T)= (1-p)^{x-1}\lambda^x e^{-\lambda t_x} * (p+(1-p)e^{-\lambda T-t_x}) = p(1-p)^{x-1}\lambda^x e^{-\lambda t_x}+(1-p)^{x}\lambda^x e^{-\lambda}L(λ,pt1,t2,...,tx,T)=(1p)x1λxeλtx(p+(1p)eλTtx)=p(1p)x1λxeλtx+(1p)xλxeλ

此时我们发现每一次交易发生的时间t是不需要的,因此我们只需要X=x(交易次数)、txt_xtx(最后一次交易时间)与T就可以了。
还有一种极端情况,就是在(0,T]之间内没有任何的交易,我们假设这用户在一开始就是活跃的,那么现在(T时刻)依然活跃的概率就是:L(λ∣X=0,T)=e−λTL(\lambda | X=0,T)=e^{-\lambda T}L(λX=0,T)=eλT(这个式子依然可以根据泊松过程发生0次事件的概率获得)
最终得到个人水平下在T时交易次数为x的似然函数

L(λ,p∣X=x,T)=δx>0p(1−p)x−1λxe−λtx+(1−p)xλxe−λL(\lambda,p | X=x,T)= \delta_{x>0} p(1-p)^{x-1}\lambda^x e^{-\lambda t_x}+(1-p)^{x}\lambda^x e^{-\lambda}L(λ,pX=x,T)=δx>0p(1p)x1λxeλtx+(1p)xλxeλ

其中δx>0\delta_{x>0}δx>0为示性函数,当x>0是为1,否则为0。

2、P(X(t)=x)P(X(t)=x)P(X(t)=x)的推导

此处的X(t)X(t)X(t)代表第t期发生交易的次数
我们将有关交易时间的变量转化为有关交易次数的变量:

X(t)>=x⇔Tx<=tX(t)>=x \Leftrightarrow T_x<=tX(t)>=xTx<=t

其中TxT_xTx表示第x次交易的时间,于是我们有:

P(X(t)=x)=P(X(t)>=x)−P(X(t)>=x+1)=P(Tx<=t)−P(Tx+1<=t)P(X(t)=x)=P(X(t)>=x)-P(X(t)>=x+1)=P(T_x<=t)-P(T_{x+1}<=t)P(X(t)=x)=P(X(t)>=x)P(X(t)>=x+1)=P(Tx<=t)P(Tx+1<=t)

而根据我们对用户失活的假设,这个式子也可以写成

P(X(t)=t)=P(第x次交易后依然活跃)P(Tx<=t且Tx+1>t)+δx>0P(第x次交易后不再活跃)P(Tx<=t)P(X(t)=t)=P(第x次交易后依然活跃)P(T_x<=t 且 T_{x+1}>t)+\delta_{x>0} P(第x次交易后不再活跃)P(T_x<=t)P(X(t)=t)=P(x次交易后依然活跃)P(Tx<=tTx+1>t)+δx>0P(x次交易后不再活跃)P(Tx<=t)

从而这个式子可以写成P(X(t)=x∣λ,p)=(1−p)xe−λt(λt)xx!+δx>0p(1−p)x−1[1−e−λtΣj=0x−1(λt)jj!]P(X(t)=x|\lambda,p)=(1-p)^x\frac{e^{-\lambda t}(\lambda t)^x}{x!} +\delta_{x>0} p(1-p)^{x-1}[1-e^{-\lambda t}\Sigma_{j=0}^{x-1}\frac{(\lambda t)^j}{j!} ]P(X(t)=xλ,p)=(1p)xx!eλt(λt)x+δx>0p(1p)x1[1eλtΣj=0x1j!(λt)j]

在原链接中,还提供了一个证明其为概率质量函数(PMF)的过程:
式子的前一项从0到无穷的累加可以写为

e−λptΣx=0∞(λ(1−p)t)xe−λ(1−p)tx!e^{-\lambda pt}\Sigma_{x=0}^\infty \frac{(\lambda(1-p)t)^x e^{-\lambda(1-p)t}}{x!}eλptΣx=0x!(λ(1p)t)xeλ(1p)t

相当于拿出一份e−λpte^{-\lambda pt}eλpt作为常数系数,而求和项就恰好是一个均值为λ(1−p)t\lambda(1-p)tλ(1p)t的泊松分布的概率质量函数,求和为1,故而整个式子的前一项等于e−λpte^{-\lambda pt}eλpt
而后一项的中括号里其实是Erlang分布的累计分布函数,故而后一项的累计求和可以写为

Σx=1∞p(1−p)x−1∫0tλxux−1e(−λu)(x−1)!du\Sigma_{x=1}^\infty p(1-p)^{x-1}\int_{0}^{t}\frac{\lambda^x u^{x-1}e^(-\lambda u)}{(x-1)!}duΣx=1p(1p)x10t(x1)!λxux1e(λu)du

=∫0tλpe−λpuΣx=1∞[λ(1−p)u]e−λ(1−p)u(x−1)!du=\int_{0}^{t}\lambda pe^{-\lambda pu}\Sigma_{x=1}^\infty \frac{[\lambda(1-p)u]e^{-\lambda(1-p)u}}{(x-1)!}du=0tλpeλpuΣx=1(x1)![λ(1p)u]eλ(1p)udu

求和项中的式子可以将x-1看做是随机变量,那么求和项就正好是一个均值为λ(1−p)u\lambda(1-p)uλ(1p)u的泊松分布的概率质量函数求和,为1,于是本式可以写为∫0tλpe−λpudu=1−e−λpt\int_{0}^{t}\lambda pe^{-\lambda pu}du=1-e^{-\lambda pt}0tλpeλpudu=1eλpt
综上,我们可以看出Σx=0∞P(X(t)=x∣λ,p)=1\Sigma_{x=0}^\infty P(X(t)=x|\lambda,p)=1Σx=0P(X(t)=xλ,p)=1,它的确是一个概率质量函数

3、P(在t时刻活跃)的推导

对于客户在t时刻还活跃的情况,有如下这么几种可能:
客户在(0,t]时刻内没有任何购买
客户在(0,t]时刻只购买了1次,而且在第1次购买后并没有“失活”
客户在(0,t]时刻只购买了2次,而且在第2次购买后并没有“失活”

从而我们能够获得如下的式子:

P(在第t时刻活跃∣λ,p)=Σj=0∞(1−p)j(λt)je−λtj!P(在第t时刻活跃|\lambda,p)=\Sigma_{j=0}^{\infty}(1-p)^j\frac{(\lambda t)^je^{-\lambda t}}{j!}P(在第t时刻活跃λ,p)=Σj=0(1p)jj!(λt)jeλt

=e−λptΣj=0∞(λ(1−p)t)je−λ(1−p)tj!=e^{-\lambda pt}\Sigma_{j=0}^{\infty}\frac{(\lambda(1-p)t)^je^{-\lambda(1-p)t}}{j!}=eλptΣj=0j!(λ(1p)t)jeλ(1p)t

上式求和部分正好是一个均值为λ(1−p)t\lambda(1-p)tλ(1p)t的破松分布的概率质量函数求和,为1。
所以P(在第t时刻活跃∣λ,p)=e−λptP(在第t时刻活跃|\lambda,p)=e^{-\lambda pt}P(在第t时刻活跃λ,p)=eλpt
我们假设用户在τ\tauτ时刻不再活跃,那么

P(τ>t∣λ,p)=e−λptP(\tau > t| \lambda,p)=e^{-\lambda pt}P(τ>tλ,p)=eλpt

根据这个定义,能够获得其累积分布函数,求导得到概率密度函数:

g(τ∣λ,p)=λpe−λpτg(\tau|\lambda,p)=\lambda pe^{-\lambda p\tau}g(τλ,p)=λpeλpτ

4、E[X(t)]的推导

和上面一样,假设用户在τ\tauτ时刻不再活跃,那么当τ>t\tau>tτ>t时,E[X(t)]就是单纯的泊松分布期望:λt\lambda tλt
而当τ<=t\tau<=tτ<=t时,在用户活跃的时间段内交易量的期望数为λτ\lambda \tauλτ

于是我们有:

E(X(t)∣λ,p)=(λt)P(τ>t)+∫0tλτg(τ∣λ,p)dτE(X(t)|\lambda,p)=(\lambda t)P(\tau>t)+\int_{0}^{t}\lambda \tau g(\tau|\lambda,p)d\tauE(X(t)λ,p)=(λt)P(τ>t)+0tλτg(τλ,p)dτ

=λte−λpt+λ2p∫0tτe−λpτdτ=\lambda te^{-\lambda pt}+\lambda^2 p\int_{0}^{t}\tau e^{-\lambda p\tau}d\tau=λteλpt+λ2p0tτeλpτdτ

=1p−1pe−λpt=\frac{1}{p}-\frac{1}{p}e^{-\lambda pt}=p1p1eλpt

对任意用户的推导

上面的所有式子都是以交易率λ\lambdaλ和用户失活率p作为条件的,但是在实际场景下这2个都是无法被观测到的。为了对任意的单个客户都能得到表达式,我们需要取λ\lambdaλ和p分布的期望。

1、似然函数推导

对于任意选择的一个客户而言,

L(r,α,a,b∣X=x,tx,T)=∫01∫0∞L(λ,p∣X=x,T)f(λ∣r,α)f(p∣a,b)dλdpL(r,\alpha,a,b|X=x,t_x,T)=\int_0^1\int_0^\infty L(\lambda,p | X=x,T) f(\lambda | r,\alpha) f(p|a,b)d\lambda dpL(r,α,a,bX=x,tx,T)=010L(λ,pX=x,T)f(λr,α)f(pa,b)dλdp

我们根据基础假设和单个客户的似然函数:

f(λ∣r,α)=αrλr−1e−λαΓ(r),λ>0f(\lambda | r,\alpha)=\frac{\alpha^r\lambda ^{r-1}e^{-\lambda\alpha}}{\Gamma(r)} , \lambda>0f(λr,α)=Γ(r)αrλr1eλα,λ>0

f(p∣a,b)=pa−1(1−p)b−1B(a,b),0<=p<=1f(p|a,b)=\frac{p^{a-1}(1-p)^{b-1}}{B(a,b)},0<=p<=1f(pa,b)=B(a,b)pa1(1p)b10<=p<=1

L(λ,p∣X=x,T)=δx>0p(1−p)x−1λxe−λtx+(1−p)xλxe−λL(\lambda,p | X=x,T)= \delta_{x>0} p(1-p)^{x-1}\lambda^x e^{-\lambda t_x}+(1-p)^{x}\lambda^x e^{-\lambda}L(λ,pX=x,T)=δx>0p(1p)x1λxeλtx+(1p)xλxeλ

原式=∫01∫0∞(δx>0p(1−p)x−1λxe−λtx+(1−p)xλxe−λ)αrλr−1e−λαΓ(r)pa−1(1−p)b−1B(a,b)dλdp\int_0^1\int_0^\infty (\delta_{x>0} p(1-p)^{x-1}\lambda^x e^{-\lambda t_x}+(1-p)^{x}\lambda^x e^{-\lambda}) \frac{\alpha^r\lambda ^{r-1}e^{-\lambda\alpha}}{\Gamma(r)} \frac{p^{a-1}(1-p)^{b-1}}{B(a,b)} d\lambda dp010(δx>0p(1p)x1λxeλtx+(1p)xλxeλ)Γ(r)αrλr1eλαB(a,b)pa1(1p)b1dλdp

注意到,

∫01pj(1−p)kpa−1(1−p)b−1B(a,b)dp=B(a+j,b+k)B(a,b)\int_0^1 p^j(1-p)^k\frac{p^{a-1}(1-p)^{b-1}}{B(a,b)}dp=\frac{B(a+j,b+k)}{B(a,b)}01pj(1p)kB(a,b)pa1(1p)b1dp=B(a,b)B(a+j,b+k)

∫0∞λje−λtαrλr−1e−λαΓ(r)dλ=Γ(r+j)αrΓ(r)(α+t)r+j\int_0^\infty \lambda^j e^{-\lambda t}\frac{\alpha^r\lambda^{r-1}e^{-\lambda\alpha}}{\Gamma(r)} d\lambda=\frac{\Gamma(r+j)\alpha^r}{\Gamma(r)(\alpha+t)^{r+j}}0λjeλtΓ(r)αrλr1eλαdλ=Γ(r)(α+t)r+jΓ(r+j)αr

这2个式子在下面会经常用到,下简称“化简式”

最终能够得到

原式=B(a,b+x)B(a,b)Γ(r+x)αrΓ(r)(α+T)r+x+δx>0B(a+1,b+x−1)B(a,b)Γ(r+x)αrΓ(r)(α+tx)r+x\frac{B(a,b+x)}{B(a,b)}\frac{\Gamma(r+x)\alpha^r}{\Gamma(r)(\alpha+T)^{r+x}} + \delta_{x>0}\frac{B(a+1,b+x-1)}{B(a,b)}\frac{\Gamma(r+x)\alpha^r}{\Gamma(r)(\alpha+t_x)^{r+x}}B(a,b)B(a,b+x)Γ(r)(α+T)r+xΓ(r+x)αr+δx>0B(a,b)B(a+1,b+x1)Γ(r)(α+tx)r+xΓ(r+x)αr

2、P(X(t)=x)的推导

对于任意选择的一个客户而言

P(X(t)=x∣r,α,a,b)=∫01∫0∞P(X(t)=x∣λ,p)f(λ∣r,α)f(p∣a,b)dλdpP(X(t)=x|r,\alpha,a,b)=\int_0^1\int_0^\infty P(X(t)=x|\lambda,p) f(\lambda | r,\alpha) f(p|a,b)d\lambda dpP(X(t)=xr,α,a,b)=010P(X(t)=xλ,p)f(λr,α)f(pa,b)dλdp

其中,

P(X(t)=x∣λ,p)=(1−p)xe−λt(λt)xx!+δx>0p(1−p)x−1[1−e−λtΣj=0x−1(λt)jj!]P(X(t)=x|\lambda,p)=(1-p)^x\frac{e^{-\lambda t}(\lambda t)^x}{x!} +\delta_{x>0} p(1-p)^{x-1}[1-e^{-\lambda t}\Sigma_{j=0}^{x-1}\frac{(\lambda t)^j}{j!} ]P(X(t)=xλ,p)=(1p)xx!eλt(λt)x+δx>0p(1p)x1[1eλtΣj=0x1j!(λt)j]

根据基础假设与化简式,可以得到

P(X(t)=x∣r,α,a,b)=B(a,b+x)B(a,b)Γ(r+x)Γ(r)x!(αα+t)r(tα+t)x+δx>0B(a+1,b+x−1)B(a,b)[1−(αα+t)r{Σj=0x−1Γ(r+j)Γ(r)j!}(tα+t)j]P(X(t)=x|r,\alpha,a,b)=\frac{B(a,b+x)}{B(a,b)}\frac{\Gamma(r+x)}{\Gamma(r)x!}(\frac{\alpha}{\alpha+t})^r(\frac{t}{\alpha+t})^x+\delta_{x>0}\frac{B(a+1,b+x-1)}{B(a,b)}[1-(\frac{\alpha}{\alpha+t})^r\left\{\Sigma_{j=0}^{x-1}\frac{\Gamma(r+j)}{\Gamma(r)j!}\right\}(\frac{t}{\alpha+t})^j]P(X(t)=xr,α,a,b)=B(a,b)B(a,b+x)Γ(r)x!Γ(r+x)(α+tα)r(α+tt)x+δx>0B(a,b)B(a+1,b+x1)[1(α+tα)r{Σj=0x1Γ(r)j!Γ(r+j)}(α+tt)j]

3、P(在t时刻依然存活)的推导

对于任意选择的一个客户而言

P(在t时刻依然存活∣r,α,a,b)=∫01∫0∞P(在第t时刻活跃∣λ,p)f(λ∣r,α)f(p∣a,b)dλdpP(在t时刻依然存活|r,\alpha,a,b)=\int_0^1\int_0^\infty P(在第t时刻活跃|\lambda,p) f(\lambda | r,\alpha) f(p|a,b) d\lambda dpP(t时刻依然存活r,α,a,b)=010P(在第t时刻活跃λ,p)f(λr,α)f(pa,b)dλdp

其中,P(在第t时刻活跃∣λ,p)=e−λptP(在第t时刻活跃|\lambda,p)=e^{-\lambda pt}P(在第t时刻活跃λ,p)=eλpt;根据我们的假设,可得:

P(在t时刻依然存活∣r,α,a,b)=∫01∫0∞e−λptαrλr−1e−λαΓ(r)pa−1(1−p)b−1B(a,b)dλdpP(在t时刻依然存活|r,\alpha,a,b)=\int_0^1\int_0^\infty e^{-\lambda pt} \frac{\alpha^r\lambda ^{r-1}e^{-\lambda\alpha}}{\Gamma(r)} \frac{p^{a-1}(1-p)^{b-1}}{B(a,b)} d\lambda dpP(t时刻依然存活r,α,a,b)=010eλptΓ(r)αrλr1eλαB(a,b)pa1(1p)b1dλdp


=αrB(a,b)∫01∫0∞λr−1e−λ(α+pt)dλΓ(r)pa−1(1−p)b−1dp=\frac{\alpha^r}{B(a,b)}\int_0^1 \frac{\int_0^\infty\lambda ^{r-1}e^{-\lambda(\alpha+pt)}d\lambda}{\Gamma(r)} p^{a-1}(1-p)^{b-1} dp=B(a,b)αr01Γ(r)0λr1eλ(α+pt)dλpa1(1p)b1dp

其中,∫0∞λr−1e−λ(α+pt)dλ=1(α+pt)r∫0∞((α+pt)λ)r−1e−λ(α+pt)d((α+pt)λ)=(α+pt)−rΓ(r)\int_0^\infty\lambda ^{r-1}e^{-\lambda(\alpha+pt)}d\lambda=\frac{1}{(\alpha+pt)^r}\int_0^\infty ((\alpha+pt)\lambda) ^{r-1}e^{-\lambda(\alpha+pt)}d((\alpha+pt)\lambda)=(\alpha+pt)^{-r}\Gamma(r)0λr1eλ(α+pt)dλ=(α+pt)r10((α+pt)λ)r1eλ(α+pt)d((α+pt)λ)=(α+pt)rΓ(r)

所以原式=αrB(a,b)∫01pa−1(1−p)b−1(α+pt)−rdp=\frac{\alpha^r}{B(a,b)}\int_0^1p^{a-1}(1-p)^{b-1}(\alpha+pt)^{-r}dp=B(a,b)αr01pa1(1p)b1(α+pt)rdp

积分换元:q=1-p,则原式=αrB(a,b)∫01qb−1(1−q)a−1(α+t−qt)−rdq=\frac{\alpha^r}{B(a,b)}\int_0^1q^{b-1}(1-q)^{a-1}(\alpha+t-qt)^{-r}dq=B(a,b)αr01qb1(1q)a1(α+tqt)rdq

=(αα+t)r1B(a,b)∫01qb−1(1−q)a−1(1−tα+tq)−rdq=(\frac{\alpha}{\alpha+t})^r\frac{1}{B(a,b)}\int_0^1q^{b-1}(1-q)^{a-1}(1-\frac{t}{\alpha+t}q)^{-r}dq=(α+tα)rB(a,b)101qb1(1q)a1(1α+ttq)rdq

而根据超几何函数的积分形式(可以看:https://zh.wikipedia.org/wiki/%E8%B6%85%E5%87%A0%E4%BD%95%E5%87%BD%E6%95%B0#%E7%A7%AF%E5%88%86%E8%A1%A8%E7%A4%BA )

2F1(a,b;c;z)=1B(b,c−b)∫01tb−1(1−t)c−b−1(1−zt)−adt,c>b_2F_1(a,b;c;z)=\frac{1}{B(b,c-b)}\int_0^1t^{b-1}(1-t)^{c-b-1}(1-zt)^{-a}dt ,\quad c>b2F1(a,b;c;z)=B(b,cb)101tb1(1t)cb1(1zt)adt,c>b

原式=P(在t时刻依然存活∣r,α,a,b)=(αα+t)r2F1(r,b;a+b;tα+t)=P(在t时刻依然存活 | r,\alpha,a,b)=\left(\frac{\alpha}{\alpha+t}\right)^r{}_2F_1\big(r,b;a+b;\frac{t}{\alpha+t}\big)=P(t时刻依然存活r,α,a,b)=(α+tα)r2F1(r,b;a+b;α+tt)

4、E[X(t)]的推导

E(X(t)∣λ,p)=1p−1pe−λptE(X(t)|\lambda,p)=\frac{1}{p}-\frac{1}{p}e^{-\lambda pt}E(X(t)λ,p)=p1p1eλpt

首先让我们求出E(X(t)∣r,α,p)E(X(t)|r,\alpha,p)E(X(t)r,α,p)。由于f(λ∣r,α)=αrλr−1e−λαΓ(r),λ>0f(\lambda | r,\alpha)=\frac{\alpha^r\lambda ^{r-1}e^{-\lambda\alpha}}{\Gamma(r)} , \lambda>0f(λr,α)=Γ(r)αrλr1eλα,λ>0,我们可以求出:

E(X(t)∣r,α,p)=1p−∫0∞1pe−λptαrλr−1e−λαΓ(r)dλE(X(t)|r,\alpha,p)=\frac{1}{p}-\int_0^\infty\frac{1}{p}e^{-\lambda pt}\frac{\alpha^r\lambda ^{r-1}e^{-\lambda\alpha}}{\Gamma(r)} d\lambdaE(X(t)r,α,p)=p10p1eλptΓ(r)αrλr1eλαdλ

右侧的积分结果在我们求P(在t时刻依然存活∣r,α,a,b)P(在t时刻依然存活|r,\alpha,a,b)P(t时刻依然存活r,α,a,b)时已经求过了,所以E(X(t)∣r,α,p)=1p−αrp(α+pt)rE(X(t)|r,\alpha,p)=\frac{1}{p}-\frac{\alpha^r}{p(\alpha+pt)^r}E(X(t)r,α,p)=p1p(α+pt)rαr
接下来我们将p的分布带入,求E(X(t)∣r,α,a,b)E(X(t)|r,\alpha,a,b)E(X(t)r,α,a,b)。其中f(p∣a,b)=pa−1(1−p)b−1B(a,b),0<=p<=1f(p|a,b)=\frac{p^{a-1}(1-p)^{b-1}}{B(a,b)},0<=p<=1f(pa,b)=B(a,b)pa1(1p)b10<=p<=1

对于被减数而言,∫011ppa−1(1−p)b−1B(a,b)dp=∫01pa−1(1−p)b−1dpB(a,b)=B(a−1)(b)B(a,b)=a+b−1a−1\int_0^1\frac1p\frac{p^{a-1}(1-p)^{b-1}}{B(a,b)}dp = \frac{\int_0^1p^{a-1}(1-p)^{b-1}dp}{B(a,b)}=\frac{B(a-1)(b)}{B(a,b)}=\frac{a+b-1}{a-1}01p1B(a,b)pa1(1p)b1dp=B(a,b)01pa1(1p)b1dp=B(a,b)B(a1)(b)=a1a+b1(贝塔函数递推公式)

对于减数而言,∫01αrp(α+pt)rpa−1(1−p)b−1B(a,b)dp=αrB(a,b)∫01pa−2(1−p)b−1(α+pt)−rdp\int_{0}^{1}\frac{\alpha^{r}}{p(\alpha+pt)^{r}}\frac{p^{a-1}(1-p)^{b-1}}{B(a,b)}dp=\frac{\alpha^r}{B(a,b)}\int_0^1p^{a-2}(1-p)^{b-1}(\alpha+pt)^{-r}dp01p(α+pt)rαrB(a,b)pa1(1p)b1dp=B(a,b)αr01pa2(1p)b1(α+pt)rdp

和上面一样,使用积分换元q=1-p,于是

减数==αrB(a,b)∫01qb−1(1−q)a−2(α+t−qt)−rdq\begin{aligned}&=\frac{\alpha^r}{B(a,b)}\int_0^1q^{b-1}(1-q)^{a-2}(\alpha+t-qt)^{-r}dq\end{aligned}=B(a,b)αr01qb1(1q)a2(α+tqt)rdq

=(αα+t)r1B(a,b)∫01qb−1(1−q)a−2(1−tα+tq)−rdq=\left(\frac{\alpha}{\alpha+t}\right)^r\frac{1}{B(a,b)}\int_0^1q^{b-1}(1-q)^{a-2}\big(1-\frac{t}{\alpha+t}q\big)^{-r}dq=(α+tα)rB(a,b)101qb1(1q)a2(1α+ttq)rdq

和上面一样,将其写作超几何函数的积分形式,得

(αα+t)rB(a−1,b)B(a,b)2F1(r,b;a+b−1;tα+t)\left(\frac{\alpha}{\alpha+t}\right)^r\frac{B(a-1,b)}{B(a,b)}_2F_1\big(r,b;a+b-1;\frac{t}{\alpha+t}\big)(α+tα)rB(a,b)B(a1,b)2F1(r,b;a+b1;α+tt)

综上,我们有

E(X(t)∣r,α,a,b)=a+b−1a−1[1−(αα+t)r2F1(r,b;a+b−1;tα+t)]E(X(t)|r,\alpha,a,b)=\frac{a+b-1}{a-1}\left[1-\left(\frac{\alpha}{\alpha+t}\right)^{r}{}_{2}F_{1}\left(r,b;a+b-1;\frac{t}{\alpha+t}\right)\right]E(X(t)r,α,a,b)=a1a+b1[1(α+tα)r2F1(r,b;a+b1;α+tt)]

参数求解

接下来我们来看一下如何求解模型中的四个参数(r,α,a,br,\alpha,a,br,α,a,b)
参考:https://brucehardie.com/notes/004/bgnbd_spreadsheet_note.pdf

根据贝塔函数的递推公式,我们有B(a+1,b+x)=aa+b+xB(a,b+x)=b+x−1a+b+xB(a+1,b+x−1)B(a+1,b+x) = \frac{a}{a+b+x}B(a,b+x) = \frac{b+x-1}{a+b+x}B(a+1,b+x-1)B(a+1,b+x)=a+b+xaB(a,b+x)=a+b+xb+x1B(a+1,b+x1)

B(a+1,b+x−1)=ab+x−1B(a,b+x)B(a+1,b+x-1) = \frac{a}{b+x-1}B(a,b+x)B(a+1,b+x1)=b+x1aB(a,b+x)

根据公式B(a,b)=Γ(a)Γ(b)Γ(a+b)B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}B(a,b)=Γ(a+b)Γ(a)Γ(b),可以得出B(a,b+x)B(a,b)=Γ(a+b)Γ(b+x)Γ(b)Γ(a+b+x)\frac{B(a,b+x)}{B(a,b)} = \frac{\Gamma(a+b)\Gamma(b+x)}{\Gamma(b)\Gamma(a+b+x)}B(a,b)B(a,b+x)=Γ(b)Γ(a+b+x)Γ(a+b)Γ(b+x)

从而似然函数能够写成L(r,α,a,b∣X=x,tx,T)=A1A2(A3+δx>0A4)L(r,\alpha,a,b|X=x,t_x,T)=A_1A_2(A_3+\delta_{x>0}A_4)L(r,α,a,bX=x,tx,T)=A1A2(A3+δx>0A4)的形式

其中,A1=Γ(r+x)αrΓ(r)A_1=\frac{\Gamma(r+x)\alpha^r}{\Gamma(r)}A1=Γ(r)Γ(r+x)αr

A2=Γ(a+b)Γ(b+x)Γ(b)Γ(a+b+x)A_2=\frac{\Gamma(a+b)\Gamma(b+x)}{\Gamma(b)\Gamma(a+b+x)}A2=Γ(b)Γ(a+b+x)Γ(a+b)Γ(b+x)

A3=(1α+T)r+xA_3=\left(\frac1{\alpha+T}\right)^{r+x}A3=(α+T1)r+x

A4=(ab+x−1)(1α+tx)r+xA_4=\Big(\frac a{b+x-1}\Big)\Big(\frac1{\alpha+t_x}\Big)^{r+x}A4=(b+x1a)(α+tx1)r+x

为了求得似然函数的对数形式,我们有

ln(A1)=ln(Γ(r+x))−ln(Γ(r))+rln(α)ln(A_1) = ln(\Gamma(r+x))-ln(\Gamma(r))+rln(\alpha)ln(A1)=ln(Γ(r+x))ln(Γ(r))+rln(α)

ln(A2)=ln(Γ(a+b))+ln(Γ(b+x))−ln(Γ(b))−ln(Γ(a+b+x))ln(A_2) = ln(\Gamma(a+b))+ln(\Gamma(b+x))-ln(\Gamma(b))-ln(\Gamma(a+b+x))ln(A2)=ln(Γ(a+b))+ln(Γ(b+x))ln(Γ(b))ln(Γ(a+b+x))

ln(A3)=−(r+x)ln(α+T)ln(A_3) = -(r+x)ln(\alpha+T)ln(A3)=(r+x)ln(α+T)

ln(A4)=ln(a)−ln(b+x−1)−(r+x)ln(α+tx)ln(A_4) = ln(a)-ln(b+x-1)-(r+x)ln(\alpha+t_x)ln(A4)=ln(a)ln(b+x1)(r+x)ln(α+tx)

最终我们有

ln[L(r,α,a,b∣X=x,tx,T)]=ln(A1)+ln(A2)+ln(A3+δx>0A4)ln[L(r,\alpha,a,b|X=x,t_x,T)]=ln(A_1)+ln(A_2)+ln(A_3+\delta_{x>0}A_4)ln[L(r,α,a,bX=x,tx,T)]=ln(A1)+ln(A2)+ln(A3+δx>0A4)

在我给出的那个页面中,作者使用了Excel作为参数求解的工具:
在这里插入图片描述
如图所示,B、C、D列3列是我们的原始数据。x代表从第0期到当前期数的总交易次数,txt_xtx代表交易期数,注意tx=0t_x=0tx=0时,x也为0,T代表真实时间与设定的第0时刻相差多少。
将4个参数放到B1到B4,然后分别在F到I列求出A1A_1A1A4A_4A4的对数值,然后在E列求出每行的积,最后在B5求出似然函数的求和。我们可以使用Excel进行规划参数求解。
在这里插入图片描述
在文件-选项中选择【加载项】,选中【规划求解加载项】,然后点积下面的【转到(G)…】
在这里插入图片描述
选中规划求解加载项,然后点击确定。
在这里插入图片描述
最后,在数据选项卡中点击右上角的规划求解按钮,注意此处的B5是似然函数的几个乘积,选择最大值。就可以使用Excel进行参数求解了。
下面我们来看一下python的lifetimes包时如何进行处理的。
这是lifetimes中BetaGeoFitter类的fit方法:

def fit(
        self, 
        frequency, 
        recency, 
        T, 
        weights=None, 
        initial_params=None, 
        verbose=False, 
        tol=1e-7, 
        index=None, 
        **kwargs
    ):
    	frequency = np.asarray(frequency).astype(int)
        recency = np.asarray(recency)
        T = np.asarray(T)
        _check_inputs(frequency, recency, T)

        if weights is None:
            weights = np.ones_like(recency, dtype=int)
        else:
            weights = np.asarray(weights)

        self._scale = _scale_time(T)
        scaled_recency = recency * self._scale
        scaled_T = T * self._scale

        log_params_, self._negative_log_likelihood_, self._hessian_ = self._fit(
            (frequency, scaled_recency, scaled_T, weights, self.penalizer_coef),
            initial_params,
            4,
            verbose,
            tol,
            **kwargs
        )

        self.params_ = pd.Series(np.exp(log_params_), index=["r", "alpha", "a", "b"])
        self.params_["alpha"] /= self._scale

        self.data = pd.DataFrame({"frequency": frequency, "recency": recency, "T": T, "weights": weights}, index=index)

        self.generate_new_data = lambda size=1: beta_geometric_nbd_model(
            T, *self._unload_params("r", "alpha", "a", "b"), size=size
        )

        self.predict = self.conditional_expected_number_of_purchases_up_to_time

        self.variance_matrix_ = self._compute_variance_matrix()
        self.standard_errors_ = self._compute_standard_errors()
        self.confidence_intervals_ = self._compute_confidence_intervals()

        return self

在所有参数中,frequency要输入顾客购买的频数,即公式的X(t)=xX(t)=xX(t)=x中的x
而recency要输入输入顾客最近一次购买的时间,即公式中的txt_xtx
最后的T即是客户的寿命(现在一直到客户第一次购买),即公式中的T
以上的三个值都可以通过lifetime中的summary_data_from_transaction_data或calibration_and_holdout_data(会留出一部分验证集)来实现。
剩下的几个参数中
weights参数可以将每个transaction加上权重(可以看做这个transaction是由几个customer购买),某个transaction的权重大于1时,表明它有多个customer,因此似然函数对数要乘以这个权重(不填默认全都是1)。
initial_params:要优化的4个参数的初始值,不填默认会给个0.1
verbose:是否要打印模型收敛的过程
tol:小于这个tolerance值之后函数终止计算
index:最终模型训练完成后会根据frequency,recency,T生成DataFrame,这个参数可以指定生成的DataFrame的索引。
在进行模型拟合之前,会根据客户寿命向量T先对T和recency进行一次缩放将其缩放到0~1之间:

def _scale_time(
    age
):
    """
    Create a scalar such that the maximum age is 1.
    """

    return 1.0 / age.max()

之后使用的是_fit方法进行的拟合:

def _fit(self, minimizing_function_args, initial_params, params_size, disp, tol=1e-7, bounds=None, **kwargs):
        # set options for minimize, if specified in kwargs will be overwritten
        minimize_options = {}
        minimize_options["disp"] = disp
        minimize_options.update(kwargs)

        current_init_params = 0.1 * np.ones(params_size) if initial_params is None else initial_params
        output = minimize(
            value_and_grad(self._negative_log_likelihood),
            jac=True,
            method=None,
            tol=tol,
            x0=current_init_params,
            args=minimizing_function_args,
            options=minimize_options,
            bounds=bounds,
        )
        if output.success:
            hessian_ = hessian(self._negative_log_likelihood)(output.x, *minimizing_function_args)
            return output.x, output.fun, hessian_
        print(output)
        raise ConvergenceError(
            dedent(
                """
            The model did not converge. Try adding a larger penalizer to see if that helps convergence.
            """
            )
        )

此处没有任何的限制条件,直接使用了BFGS方法进行优化。BFGS算法我个人理解就是使用一个正定矩阵B_k来近似Heissen矩阵以此弥补牛顿法中需要求函数二阶导的缺陷。(优化算法笔者不熟,有误的话评论区还请多多指正)。

模型预测

1、在T时刻用户依然活跃的概率

也就是P(在T时刻依然活跃∣X=x,tx,T)P(在T时刻依然活跃|X=x,t_x,T)P(T时刻依然活跃X=x,tx,T)
针对于以此都没有购买的客户(x=0),我们假定P(在T时刻依然活跃T∣X=0,T,r,α,a,b)=1P(\text{在T时刻依然活跃}T\mid X=0,T,r,\alpha,a,b)=1P(T时刻依然活跃TX=0,T,r,α,a,b)=1
而对于购买次数为x的客户,在T时刻活跃的概率为:

P(在T时刻依然活跃∣X=x,tx,T,λ,p)=(1−p)e−λ(T−tx)p+(1−p)e−λ(T−tx)P(在T时刻依然活跃\mid X=x,t_x,T,\lambda,p)=\frac{(1-p)e^{-\lambda(T-t_x)}}{p+(1-p)e^{-\lambda(T-t_x)}}P(T时刻依然活跃X=x,tx,T,λ,p)=p+(1p)eλ(Ttx)(1p)eλ(Ttx)

其中,p为用户失活的情形,而(1−p)e−λ(T−tx)(1-p)e^{-\lambda(T-t_x)}(1p)eλ(Ttx)则代表了用户在x次购买之后,在最后一次交易直到T时间没有任何交易,但是依然活跃的概率。

分子分母同时乘以(1−p)x−1λxe−λtx(1-p)^{x-1}\lambda^xe^{-\lambda t_x}(1p)x1λxeλtx

其中分母为:

L(λ,p∣X=x,tx,T)=δx>0p(1−p)x−1λxe−λtx+(1−p)xλxe−λL(\lambda,p | X=x,t_x,T)= \delta_{x>0} p(1-p)^{x-1}\lambda^x e^{-\lambda t_x}+(1-p)^{x}\lambda^x e^{-\lambda}L(λ,pX=x,tx,T)=δx>0p(1p)x1λxeλtx+(1p)xλxeλ(当x=0时,原式=1,和此处的δx>0\delta_{x>0}δx>0没有冲突)

而分子=(1−p)xλxe−λT(1-p)^{x}\lambda^xe^{-\lambda T}(1p)xλxeλT

由于式子中含有不可测的pppλ\lambdaλ,我们还需要进行计算。和前面的对任意的单个客户求表达式的思路一样,我们求:

P(在T时刻依然活跃∣X=x,tx,T,r,α,a,b)=∫01∫0∞P(在T时刻依然活跃∣X=x,tx,T,λ,p)f(λ,p∣r,α,a,b,X=x,tx,T)dλdpP(在T时刻依然活跃\mid X=x,t_x,T,r,\alpha,a,b)=\int_0^1\int_0^\infty P(在T时刻依然活跃\mid X=x,t_x,T,\lambda,p) f(\lambda,p|r,\alpha,a,b,X=x,t_x,T)d\lambda dpP(T时刻依然活跃X=x,tx,T,r,α,a,b)=010P(T时刻依然活跃X=x,tx,T,λ,p)f(λ,pr,α,a,b,X=x,tx,T)dλdp

而根据贝叶斯公式,我们有:
f(λ,p∣r,α,a,b,X=x,tx,T)=L(λ,p∣X=x,tx,T)f(λ∣r,α)f(p∣a,b)L(r,α,a,b∣X=x,tx,T)f(\lambda,p | r,\alpha,a,b,X=x,t_x,T)=\frac{L(\lambda,p | X=x,t_x,T)f(\lambda | r,\alpha)f(p | a,b)}{L(r,\alpha,a,b | X=x,t_x,T)}f(λ,pr,α,a,b,X=x,tx,T)=L(r,α,a,bX=x,tx,T)L(λ,pX=x,tx,T)f(λr,α)f(pa,b)

我们将贝叶斯公式得出的式子和我们一开始求的L(λ,p∣X=x,tx,T)L(\lambda,p | X=x,t_x,T)L(λ,pX=x,tx,T)带入到中间的P(在T时刻依然活跃∣X=x,tx,T,r,α,a,b)P(在T时刻依然活跃\mid X=x,t_x,T,r,\alpha,a,b)P(T时刻依然活跃X=x,tx,T,r,α,a,b)中,可得

P(在T时刻依然活跃∣X=x,tx,T,r,α,a,b)=1L(r,α,a,b∣X=x,tx,T)∫01∫0∞(1−p)xλxe−λTf(λ∣r,α)f(p∣a,b)dλdpP(在T时刻依然活跃\mid X=x,t_x,T,r,\alpha,a,b)=\frac{1}{L(r,\alpha,a,b | X=x,t_x,T)}\int_0^1\int_0^\infty (1-p)^x\lambda^xe^{-\lambda T}f(\lambda | r,\alpha)f(p | a,b)d\lambda dpP(T时刻依然活跃X=x,tx,T,r,α,a,b)=L(r,α,a,bX=x,tx,T)1010(1p)xλxeλTf(λr,α)f(pa,b)dλdp

式中的积分部分可以写为∫01(1−p)xf(p∣a,b)dp∫0∞λxe−λTf(λ∣r,α)dλ\int_0^1 (1-p)^xf(p | a,b) dp\int_0^\infty\lambda^xe^{-\lambda T}f(\lambda | r,\alpha)d\lambda01(1p)xf(pa,b)dp0λxeλTf(λr,α)dλ
这2个形式正好对应了我们在求单个客户表达式时使用的“化简式”,于是我们有:

∫01∫0∞(1−p)xλxe−λTf(λ∣r,α)f(p∣a,b)dλdp=B(a,b+x)B(a,b)=Γ(r+x)αrΓ(r)(α+T)r+x\int_0^1\int_0^\infty (1-p)^x\lambda^xe^{-\lambda T}f(\lambda | r,\alpha)f(p | a,b)d\lambda dp=\frac{B(a,b+x)}{B(a,b)}=\frac{\Gamma(r+x)\alpha^r}{\Gamma(r)(\alpha+T)^{r+x}}010(1p)xλxeλTf(λr,α)f(pa,b)dλdp=B(a,b)B(a,b+x)=Γ(r)(α+T)r+xΓ(r+x)αr

最终再把L(r,α,a,b∣X=x,tx,T)=B(a,b+x)B(a,b)Γ(r+x)αrΓ(r)(α+T)r+x+δx>0B(a+1,b+x−1)B(a,b)Γ(r+x)αrΓ(r)(α+tx)r+xL(r,\alpha,a,b|X=x,t_x,T)=\frac{B(a,b+x)}{B(a,b)}\frac{\Gamma(r+x)\alpha^r}{\Gamma(r)(\alpha+T)^{r+x}} + \delta_{x>0}\frac{B(a+1,b+x-1)}{B(a,b)}\frac{\Gamma(r+x)\alpha^r}{\Gamma(r)(\alpha+t_x)^{r+x}}L(r,α,a,bX=x,tx,T)=B(a,b)B(a,b+x)Γ(r)(α+T)r+xΓ(r+x)αr+δx>0B(a,b)B(a+1,b+x1)Γ(r)(α+tx)r+xΓ(r+x)αr带入,可以得:

P(在T时刻依然活跃∣X=x,tx,T,r,α,a,b)=11+δx>0ab+x−1(α+Tα+tx)r+xP(在T时刻依然活跃\mid X=x,t_x,T,r,\alpha,a,b)=\frac{1}{1+\delta_{x>0}\frac{a}{b+x-1}\left(\frac{\alpha+T}{\alpha+t_{x}}\right)^{r+x}}P(T时刻依然活跃X=x,tx,T,r,α,a,b)=1+δx>0b+x1a(α+txα+T)r+x1

2、预测E(Y(t)∣X=x,tx,T)E(Y(t)|X=x,t_x,T)E(Y(t)X=x,tx,T)

式中的Y(t)代表在时间段(T,T+t]的交易数

在最初,我们已经推导了E(Y(t)∣λ,p)=1p−1pe−λptE(Y(t)|\lambda,p)=\frac{1}{p}-\frac{1}{p}e^{-\lambda pt}E(Y(t)λ,p)=p1p1eλpt

将这个式子乘以P(在T时刻依然活跃∣X=x,tx,T)=(1−p)xλxe−λTL(λ,p∣X=x,tx,T)P(在T时刻依然活跃\mid X=x,t_x,T)=\frac{(1-p)^{x}\lambda^xe^{-\lambda T}}{L(\lambda,p | X=x,t_x,T)}P(T时刻依然活跃X=x,tx,T)=L(λ,pX=x,tx,T)(1p)xλxeλT

E(Y(t)∣X=x,tx,T,λ,p)E(Y(t)|X=x,t_x,T,\lambda,p)E(Y(t)X=x,tx,T,λ,p)

=(1−p)xλxe−λT(1p−1pe−λpt)L(λ,p∣X=x,tx,T)=\frac{(1-p)^{x}\lambda^{x}e^{-\lambda T}\left(\frac{1}{p}-\frac{1}{p}e^{-\lambda pt}\right)}{L(\lambda,p\mid X=x,t_{x},T)}=L(λ,pX=x,tx,T)(1p)xλxeλT(p1p1eλpt)

=p−1(1−p)xλxe−λT−p−1(1−p)xλxe−λ(T+pt)L(λ,p∣X=x,tx,T)=\frac{p^{-1}(1-p)^{x}\lambda^{x}e^{-\lambda T}-p^{-1}(1-p)^{x}\lambda^{x}e^{-\lambda(T+pt)}}{L(\lambda,p | X=x,t_{x},T)}=L(λ,pX=x,tx,T)p1(1p)xλxeλTp1(1p)xλxeλ(T+pt)

和上面一样,我们求E(Y(t)∣X=x,tx,T,r,α,a,b)=∫01∫0∞E(Y(t)∣X=x,tx,T,λ,p)f(λ,p∣r,α,a,b,X=x,tx,T)dλdpE(Y(t)|X=x,t_x,T,r,\alpha,a,b)=\int_0^1\int_0^\infty E(Y(t)|X=x,t_x,T,\lambda,p) f(\lambda,p|r,\alpha,a,b,X=x,t_x,T)d\lambda dpE(Y(t)X=x,tx,T,r,α,a,b)=010E(Y(t)X=x,tx,T,λ,p)f(λ,pr,α,a,b,X=x,tx,T)dλdp

=A−BL(r,α,a,b∣X=x,tx,T)=\frac{\mathrm{A}-\mathrm{B}}{L(r,\alpha,a,b\mid X=x,t_x,T)}=L(r,α,a,bX=x,tx,T)AB

其中,A=∫01∫0∞p−1(1−p)xλxe−λTf(λ∣r,α)f(p∣a,b)dλdpA=\int_0^1\int_0^\infty p^{-1}(1-p)^x\lambda^xe^{-\lambda T}f(\lambda | r,\alpha)f(p | a,b)d\lambda dpA=010p1(1p)xλxeλTf(λr,α)f(pa,b)dλdp

(根据“化简式”)=B(a−1,b+x)B(a,b)=Γ(r+x)αrΓ(r)(α+T)r+x=\frac{B(a-1,b+x)}{B(a,b)}=\frac{\Gamma(r+x)\alpha^r}{\Gamma(r)(\alpha+T)^{r+x}}=B(a,b)B(a1,b+x)=Γ(r)(α+T)r+xΓ(r+x)αr

而B=∫01∫0∞p−1(1−p)xλxe−λ(T+pt)f(λ∣r,α)f(p∣a,b)dλdp=\int_{0}^{1}\int_{0}^{\infty}p^{-1}(1-p)^{x}\lambda^{x}e^{-\lambda(T+pt)}f(\lambda | r,\alpha)f(p | a,b)d\lambda dp=010p1(1p)xλxeλ(T+pt)f(λr,α)f(pa,b)dλdp

=∫01pa−2(1−p)b+x−1B(a,b){∫0∞αrλr+x−1e−λ(α+T+pt)Γ(r)dλ}dp=\int_0^1\frac{p^{a-2}(1-p)^{b+x-1}}{B(a,b)}\left\{\int_0^\infty\frac{\alpha^r\lambda^{r+x-1}e^{-\lambda(\alpha+T+pt)}}{\Gamma(r)}d\lambda\right\}dp=01B(a,b)pa2(1p)b+x1{0Γ(r)αrλr+x1eλ(α+T+pt)dλ}dp

其中关于λ\lambdaλ的积分可以写作

∫0∞αrλr+x−1e−λ(α+T+pt)Γ(r)dλ=αrΓ(r)∫0∞ λr+x−1e−λ(α+T+pt)dλ=αrΓ(r)(α+T+pt)−(r+x)Γ(r+x)\int_0^\infty\frac{\alpha^r\lambda^{r+x-1}e^{-\lambda(\alpha+T+pt)}}{\Gamma(r)}d\lambda=\frac{\alpha^r}{\Gamma(r)}\int_0^\infty\ \lambda^{r+x-1}e^{-\lambda(\alpha+T+pt)}d\lambda=\frac{\alpha^r}{\Gamma(r)}(\alpha+T+pt)^{-(r+x)}\Gamma(r+x)0Γ(r)αrλr+x1eλ(α+T+pt)dλ=Γ(r)αr0 λr+x1eλ(α+T+pt)dλ=Γ(r)αr(α+T+pt)(r+x)Γ(r+x)

从而B=Γ(r+x)αrΓ(r)B(a,b)∫01pa−2(1−p)b+x−1(α+T+pt)−(r+x)dpB=\frac{\Gamma(r+x)\alpha^r}{\Gamma(r)B(a,b)}\int_0^1p^{a-2}(1-p)^{b+x-1}(\alpha+T+pt)^{-(r+x)}dpB=Γ(r)B(a,b)Γ(r+x)αr01pa2(1p)b+x1(α+T+pt)(r+x)dp

和上面一样,使用q=1−pq=1-pq=1p进行还元,得

Γ(r+x)αrΓ(r)B(a,b)∫01qb+x−1(1−q)a−2(α+T+t−qt)−(r+x)dq\frac{\Gamma(r+x)\alpha^r}{\Gamma(r)B(a,b)}\int_0^1q^{b+x-1}(1-q)^{a-2}(\alpha+T+t-qt)^{-(r+x)}dqΓ(r)B(a,b)Γ(r+x)αr01qb+x1(1q)a2(α+T+tqt)(r+x)dq

=Γ(r+x)αrΓ(r)B(a,b)(α+T+t)r+x∫01qb+x−1(1−q)a−2(1−tα+T+tq)−(r+x)dq=\frac{\Gamma(r+x)\alpha^r}{\Gamma(r)B(a,b)(\alpha+T+t)^{r+x}}\int_0^1q^{b+x-1}(1-q)^{a-2}\big(1-\frac t{\alpha+T+t}q\big)^{-(r+x)}dq=Γ(r)B(a,b)(α+T+t)r+xΓ(r+x)αr01qb+x1(1q)a2(1α+T+ttq)(r+x)dq

积分部分可以写作超几何函数的积分形式,故而

B=B(a−1,b+x)B(a,b)Γ(r+x)αrΓ(r)(α+T+t)r+x2F1(r+x,b+x;a+b+x−1;tα+T+t)B=\frac{B(a-1,b+x)}{B(a,b)}\frac{\Gamma(r+x)\alpha^{r}}{\Gamma(r)(\alpha+T+t)^{r+x}}_2F_1\big(r+x,b+x;a+b+x-1;\frac{t}{\alpha+T+t}\big)B=B(a,b)B(a1,b+x)Γ(r)(α+T+t)r+xΓ(r+x)αr2F1(r+x,b+x;a+b+x1;α+T+tt)

综上,我们可以求得:

E(Y(t)∣X=x,tx,T,r,α,a,b)=a+b+x−1a−1[1−(α+Tα+T+t)r+x2F1(r+x,b+x;a+b+x−1;tα+T+t)]1+δx>0ab+x−1(α+Tα+tx)r+x\begin{aligned}&E(Y(t)\mid X=x,t_x,T,r,\alpha,a,b)=\\&\frac{\frac{a+b+x-1}{a-1}\left[1-\left(\frac{\alpha+T}{\alpha+T+t}\right)^{r+x} 2F_1\big(r+x,b+x;a+b+x-1;\frac{t}{\alpha+T+t}\big)\right]}{1+\delta_{x>0}\frac{a}{b+x-1}\left(\frac{\alpha+T}{\alpha+t_x}\right)^{r+x}}\end{aligned}E(Y(t)X=x,tx,T,r,α,a,b)=1+δx>0b+x1a(α+txα+T)r+xa1a+b+x1[1(α+T+tα+T)r+x2F1(r+x,b+x;a+b+x1;α+T+tt)]

3、预测P(Y(t)=y∣X=x,yx,T)P(Y(t)=y|X=x,y_x,T)P(Y(t)=yX=x,yx,T)

也就是根据一个客户的历史购买情况预测出这个客户在(T,T+t]时间段内的购买可能。
首先我们之前已经推导过了当用户在T时刻依然活跃的情形:

P(X(t)=x∣λ,p)=(1−p)xe−λt(λt)xx!+δx>0p(1−p)x−1[1−e−λtΣj=0x−1(λt)jj!]P(X(t)=x|\lambda,p)=(1-p)^x\frac{e^{-\lambda t}(\lambda t)^x}{x!} +\delta_{x>0} p(1-p)^{x-1}[1-e^{-\lambda t}\Sigma_{j=0}^{x-1}\frac{(\lambda t)^j}{j!}]P(X(t)=xλ,p)=(1p)xx!eλt(λt)x+δx>0p(1p)x1[1eλtΣj=0x1j!(λt)j]
而对于用户在T时刻不再活跃的情情形,我们有

P(Y(t)=y∣λ,p)={1如果y=00其他情况P(Y(t)=y | \lambda,p)=\begin{cases}1&\text{如果}y=0\\0&\text{其他情况}\end{cases}P(Y(t)=yλ,p)={10如果y=0其他情况

我们分别将P(在T时刻依然活跃∣X=x,tx,T,λ,p)=(1−p)xλxe−λTL(λ,p∣X=x,tx,T)P(在T时刻依然活跃|X=x,t_x,T,\lambda,p)=\frac{(1-p)^{x}\lambda^xe^{-\lambda T}}{L(\lambda,p | X=x,t_x,T)}P(T时刻依然活跃X=x,tx,T,λ,p)=L(λ,pX=x,tx,T)(1p)xλxeλT


P(在tx时刻活跃,但是到了第x次交易之后立刻失活)=p(1−p)x−1λxe−λtxL(λ,p∣X=x,tx,T)P(在t_x时刻活跃,但是到了第x次交易之后立刻失活)=\frac{p(1-p)^{x-1}\lambda^xe^{-\lambda t_x}}{L(\lambda,p | X=x,t_x,T)}P(tx时刻活跃,但是到了第x次交易之后立刻失活)=L(λ,pX=x,tx,T)p(1p)x1λxeλtx

作为系数,乘以这2个情形下的式子,获得统一的表达式:

P(Y(t)=y∣X=x,tx,T,λ,p)={δx>0,y=0p(1−p)x−1λxe−λtx+(1−p)x+yλx+ytye−λ(T+t)y!+δy>0p(1−p)x+y−1[λxe−λT−e−λ(T+t)∑j=0y−1λx+jtjj!]}/L(λ,p∣X=x,tx,T)P(Y(t)=y | X=x,t_x,T,\lambda,p)=\left\{\delta_{x>0,y=0} p(1-p)^{x-1}\lambda^xe^{-\lambda t_x}+(1-p)^{x+y}\frac{\lambda^{x+y}t^ye^{-\lambda(T+t)}}{y!}+\delta_{y>0}p(1-p)^{x+y-1}\left[\lambda^xe^{-\lambda T}-e^{-\lambda(T+t)}\sum_{j=0}^{y-1}\frac{\lambda^{x+j}t^j}{j!}\right]\right\}\Big/L(\lambda,p | X=x,t_x,T)P(Y(t)=yX=x,tx,T,λ,p)={δx>0,y=0p(1p)x1λxeλtx+(1p)x+yy!λx+ytyeλ(T+t)+δy>0p(1p)x+y1[λxeλTeλ(T+t)j=0y1j!λx+jtj]}/L(λ,pX=x,tx,T)

和上面一样,我们需将λ\lambdaλ和p转换为我们求出的那4个参数:

P(Y(t)=y∣X=x,tx,T,r,α,a,b)=∫01∫0∞{P(Y(t)=y∣X=x,tx,T,λ,p)f(λ,p∣r,α,a,b,X=x,tx,T)}dλdp\begin{aligned}P(Y(t)=y | X&=x,t_x,T,r,\alpha,a,b)\\&=\int_0^1\int_0^\infty\Big\{P(Y(t)=y | X=x,t_x,T,\lambda,p) f(\lambda,p | r,\alpha,a,b,X=x,t_x,T)\Big\} d\lambda dp\end{aligned}P(Y(t)=yX=x,tx,T,r,α,a,b)=010{P(Y(t)=yX=x,tx,T,λ,p)f(λ,pr,α,a,b,X=x,tx,T)}dλdp

这个$ f(\lambda,p | r,\alpha,a,b,X=x,t_x,T)$我们上文已经求过了,现在能够得到:

P(Y(t)=y∣X=x,tx,T,r,α,a,b)=δx>0,y=0A+B+δy>0CL(r,α,a,b∣X=x,tx,T)P(Y(t)=y | X=x,t_x,T,r,\alpha,a,b)=\frac{\delta_{x>0,y=0}\text{A}+\text{B}+\delta_{y>0}\text{C}}{L(r,\alpha,a,b | X=x,t_x,T)}P(Y(t)=yX=x,tx,T,r,α,a,b)=L(r,α,a,bX=x,tx,T)δx>0,y=0A+B+δy>0C

其中,A=B(a+1,b+x−1)B(a,b)Γ(r+x)Γ(r)αr(α+tx)r+xA=\frac{B(a+1,b+x-1)}{B(a,b)}\frac{\Gamma(r+x)}{\Gamma(r)}\frac{\alpha^r}{(\alpha+t_x)^{r+x}}A=B(a,b)B(a+1,b+x1)Γ(r)Γ(r+x)(α+tx)r+xαr
,

B=B(a,b+x+y)B(a,b)Γ(r+x+y)Γ(r)y!αrty(α+T+t)r+x+yB=\frac{B(a,b+x+y)}{B(a,b)}\frac{\Gamma(r+x+y)}{\Gamma(r)y!}\frac{\alpha^rt^y}{(\alpha+T+t)^{r+x+y}}B=B(a,b)B(a,b+x+y)Γ(r)y!Γ(r+x+y)(α+T+t)r+x+yαrty

C=B(a+1,b+x+y−1)B(a,b){Γ(r+x)Γ(r)αr(α+T)r+x−∑j=0y−1Γ(r+x+j)Γ(r)j!αrtj(α+T+t)r+x+j}C=\frac{B(a+1,b+x+y-1)}{B(a,b)}\left\{\begin{aligned}\frac{\Gamma(r+x)}{\Gamma(r)}\frac{\alpha^r}{(\alpha+T)^{r+x}}-\sum_{j=0}^{y-1}\frac{\Gamma(r+x+j)}{\Gamma(r)j!}\frac{\alpha^rt^j}{(\alpha+T+t)^{r+x+j}}\end{aligned}\right\}C=B(a,b)B(a+1,b+x+y1)Γ(r)Γ(r+x)(α+T)r+xαrj=0y1Γ(r)j!Γ(r+x+j)(α+T+t)r+x+jαrtj

lifetimes中的函数

在lifetimes包中,使用conditional_expected_number_of_purchases_up_to_time函数来预测t期之后的交易数,也就是我们求的E(Y(t)∣X=x,tx,T)E(Y(t)|X=x,t_x,T)E(Y(t)X=x,tx,T)

def conditional_expected_number_of_purchases_up_to_time(
    self, 
    t, 
    frequency, 
    recency, 
    T
):
    """
    Conditional expected number of purchases up to time.

    Calculate the expected number of repeat purchases up to time t for a
    randomly chosen individual from the population, given they have
    purchase history (frequency, recency, T).

    This function uses equation (10) from [2]_.

    Parameters
    ----------
    t: array_like
        times to calculate the expectation for.
    frequency: array_like
        historical frequency of customer.
    recency: array_like
        historical recency of customer.
    T: array_like
        age of the customer.

    Returns
    -------
    array_like

    References
    ----------
    .. [2] Fader, Peter S., Bruce G.S. Hardie, and Ka Lok Lee (2005a),
    "Counting Your Customers the Easy Way: An Alternative to the
    Pareto/NBD Model," Marketing Science, 24 (2), 275-84.
    """

    x = frequency
    r, alpha, a, b = self._unload_params("r", "alpha", "a", "b")

    _a = r + x
    _b = b + x
    _c = a + b + x - 1
    _z = t / (alpha + T + t)
    ln_hyp_term = np.log(hyp2f1(_a, _b, _c, _z))

    # if the value is inf, we are using a different but equivalent
    # formula to compute the function evaluation.
    ln_hyp_term_alt = np.log(hyp2f1(_c - _a, _c - _b, _c, _z)) + (_c - _a - _b) * np.log(1 - _z)
    ln_hyp_term = np.where(np.isinf(ln_hyp_term), ln_hyp_term_alt, ln_hyp_term)
    first_term = (a + b + x - 1) / (a - 1)
    second_term = 1 - np.exp(ln_hyp_term + (r + x) * np.log((alpha + T) / (alpha + t + T)))

    numerator = first_term * second_term
    denominator = 1 + (x > 0) * (a / (b + x - 1)) * ((alpha + T) / (alpha + recency)) ** (r + x)

    return numerator / denominator

它的式子基本上和我们求的结果一样,唯一的不同点在于分子部分带有超越函数的那部分会先进行对数化处理,如果对数化的结果为inf,则使用不同但等效的公式(ln_hyp_term_alt)来计算函数的求值。

接下来是计算在T时刻依然活跃的函数。

def conditional_probability_alive(
    self, 
    frequency, 
    recency, 
    T
):
    """
    Compute conditional probability alive.

    Compute the probability that a customer with history
    (frequency, recency, T) is currently alive.

    From http://www.brucehardie.com/notes/021/palive_for_BGNBD.pdf

    Parameters
    ----------
    frequency: array or scalar
        historical frequency of customer.
    recency: array or scalar
        historical recency of customer.
    T: array or scalar
        age of the customer.

    Returns
    -------
    array
        value representing a probability
    """

    r, alpha, a, b = self._unload_params("r", "alpha", "a", "b")

    log_div = (r + frequency) * np.log((alpha + T) / (alpha + recency)) + np.log(
        a / (b + np.maximum(frequency, 1) - 1)
    )

    return np.atleast_1d(np.where(frequency == 0, 1.0, expit(-log_div)))

def conditional_probability_alive_matrix(
    self, 
    max_frequency=None, 
    max_recency=None
):
    """
    Compute the probability alive matrix.

    Uses the ``conditional_probability_alive()`` method to get calculate the matrix.

    Parameters
    ----------
    max_frequency: float, optional
        the maximum frequency to plot. Default is max observed frequency.
    max_recency: float, optional
        the maximum recency to plot. This also determines the age of the
        customer. Default to max observed age.

    Returns
    -------
    matrix:
        A matrix of the form [t_x: historical recency, x: historical frequency]
    """

    max_frequency = max_frequency or int(self.data["frequency"].max())
    max_recency = max_recency or int(self.data["T"].max())

    return np.fromfunction(
        self.conditional_probability_alive, (max_frequency + 1, max_recency + 1), T=max_recency
    ).T

公式和我们推出来的一样,只是计算的是式子的对数形式,最后再指数展开。
除此以外,还会有计算t期内交易数期望值的函数,也就是我们计算的E(X(t)∣r,α,a,b)E(X(t)|r,\alpha,a,b)E(X(t)r,α,a,b)

def expected_number_of_purchases_up_to_time(
    self, 
    t
):
    """
    Calculate the expected number of repeat purchases up to time t.

    Calculate repeat purchases for a randomly chosen individual from the
    population.

    Equivalent to equation (9) of [2]_.

    Parameters
    ----------
    t: array_like
        times to calculate the expection for

    Returns
    -------
    array_like

    References
    ----------
    .. [2] Fader, Peter S., Bruce G.S. Hardie, and Ka Lok Lee (2005a),
    "Counting Your Customers the Easy Way: An Alternative to the
    Pareto/NBD Model," Marketing Science, 24 (2), 275-84.
    """

    r, alpha, a, b = self._unload_params("r", "alpha", "a", "b")
    hyp = hyp2f1(r, b, a + b - 1, t / (alpha + t))

    return (a + b - 1) / (a - 1) * (1 - hyp * (alpha / (alpha + t)) ** r)

公式和我们推导出来的结果一样。

最后,还有一个可以预测未来t期,交易数为n概率的函数,也就是我们的P(Y(t)=y∣X=x,yx,T)P(Y(t)=y|X=x,y_x,T)P(Y(t)=yX=x,yx,T)

def probability_of_n_purchases_up_to_time(
    self, 
    t, 
    n
):
    r"""
    Compute the probability of n purchases.

     .. math::  P( N(t) = n | \text{model} )

    where N(t) is the number of repeat purchases a customer makes in t
    units of time.

    Comes from equation (8) of [2]_.

    Parameters
    ----------
    t: float
        number units of time
    n: int
        number of purchases

    Returns
    -------
    float:
        Probability to have n purchases up to t units of time

    References
    ----------
    .. [2] Fader, Peter S., Bruce G.S. Hardie, and Ka Lok Lee (2005a),
    "Counting Your Customers the Easy Way: An Alternative to the
    Pareto/NBD Model," Marketing Science, 24 (2), 275-84.
    """

    r, alpha, a, b = self._unload_params("r", "alpha", "a", "b")

    first_term = (
        beta(a, b + n)
        / beta(a, b)
        * gamma(r + n)
        / gamma(r)
        / gamma(n + 1)
        * (alpha / (alpha + t)) ** r
        * (t / (alpha + t)) ** n
    )

    if n > 0:
        j = np.arange(0, n)
        finite_sum = (gamma(r + j) / gamma(r) / gamma(j + 1) * (t / (alpha + t)) ** j).sum()
        second_term = beta(a + 1, b + n - 1) / beta(a, b) * (1 - (alpha / (alpha + t)) ** r * finite_sum)
    else:
        second_term = 0

    return first_term + second_term
    

形式上的和我们推出来的公式不大一样,但是结果应该是一样的,注意此处的参数n应该代表的是我们推导中的x+y

总结

本文跟随了https://www.brucehardie.com/notes/039/bgnbd_derivation__2019-11-06.pdf 中的过程一步步推导了BG/NBD的数学公式,并介绍了使用Excel与lifetimes包中的求解以及预测方法。笔者本人也跟随着扩展了贝塔/伽马/泊松/指数分布的知识。之后有机会会做一个使用BG/NBD的CLV分析案例。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值