模型推导：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(t_j|t_{j-1};\lambda )=\lambda e^{-\lambda (t_j-t_{j-1})}$ , t_j>t_{j-1}>=0

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

$f(\lambda |r,\alpha )=\frac{{\alpha }^r{\lambda }^{r-1}{e}^{-\lambda \alpha }}{\Gamma (r)},\lambda >0$

其中$\Gamma (r)={\int }_0^{+\infty }{t}^{r-1}{e}^{-t}dt$，是gamma函数

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

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

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

$f(p|a,b)=\frac{p^{a-1}(1-p{)}^{b-1}}{B(a,b)}，0<=p<=1$

其中，B(a,b)是贝塔函数，$B(a,b)=\frac{\Gamma (a)\Gamma (b)}{\Gamma (a+b)}={\int }_0^1{x}^{a-1}{\left(1-x\right)}^{b-1}dx$

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

对于某个特定顾客的推导（已知$\lambda$和p）

1、似然函数的推导

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

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

$L(\lambda ,p|t_1,t_2,...,t_x,T)=(1-p{)}^{x-1}{\lambda }^x{e}^{-\lambda t_x}\ast (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 }$

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

$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 }$

其中${\delta }_{x>0}$为示性函数，当x>0是为1，否则为0。

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

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

$X(t)>=x\iff T_x<=t$

其中$T_x$表示第x次交易的时间，于是我们有:

$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)=t)=P(第x次交易后依然活跃)P(T_x<=t且T_{x+1}>t)+{\delta }_{x>0}P(第x次交易后不再活跃)P(T_x<=t)$

从而这个式子可以写成$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!}]$

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

$e^{-\lambda pt}{\Sigma }_{x=0}^{\infty }\frac{(\lambda (1-p)t{)}^x{e}^{-\lambda (1-p)t}}{x!}$

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

${\Sigma }_{x=1}^{\infty }p(1-p{)}^{x-1}{\int }_0^t\frac{{\lambda }^x{u}^{x-1}{e}^{(}-\lambda u)}{(x-1)!}du$

$={\int }_0^t\lambda p{e}^{-\lambda pu}{\Sigma }_{x=1}^{\infty }\frac{[\lambda (1-p)u]e^{-\lambda (1-p)u}}{(x-1)!}du$

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

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

对于客户在t时刻还活跃的情况，有如下这么几种可能：
客户在(0,t]时刻内没有任何购买
客户在(0,t]时刻只购买了1次，而且在第1次购买后并没有“失活”
客户在(0,t]时刻只购买了2次，而且在第2次购买后并没有“失活”
…
从而我们能够获得如下的式子：

$P(在第t时刻活跃|\lambda ,p)={\Sigma }_{j=0}^{\infty }(1-p{)}^j\frac{(\lambda t{)}^j{e}^{-\lambda t}}{j!}$

$=e^{-\lambda pt}{\Sigma }_{j=0}^{\infty }\frac{(\lambda (1-p)t{)}^j{e}^{-\lambda (1-p)t}}{j!}$

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

$P(\tau >t|\lambda ,p)=e^{-\lambda pt}$

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

$g(\tau |\lambda ,p)=\lambda p{e}^{-\lambda p\tau }$

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

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

于是我们有：

$E(X(t)|\lambda ,p)=(\lambda t)P(\tau >t)+{\int }_0^t\lambda \tau g(\tau |\lambda ,p)d\tau$

$=\lambda t{e}^{-\lambda pt}+{\lambda }^{2}p{\int }_0^t\tau e^{-\lambda p\tau }d\tau$

$=\frac{1}{p}-\frac{1}{p}{e}^{-\lambda pt}$

对任意用户的推导

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

1、似然函数推导

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

$L(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 dp$

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

$f(\lambda |r,\alpha )=\frac{{\alpha }^r{\lambda }^{r-1}{e}^{-\lambda \alpha }}{\Gamma (r)},\lambda >0$

$f(p|a,b)=\frac{p^{a-1}(1-p{)}^{b-1}}{B(a,b)}，0<=p<=1$

$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 }$

原式=${\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 dp$

注意到，

${\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)}$

${\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}}$

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

最终能够得到

原式=$\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}}$

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

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

$P(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 dp$

其中，

$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|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]$

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

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

$P(在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 dp$

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

$P(在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 dp$


$=\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$

其中，${\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)$

所以原式$=\frac{{\alpha }^r}{B(a,b)}{\int }_0^1{p}^{a-1}(1-p{)}^{b-1}(\alpha +pt{)}^{-r}dp$

积分换元:q=1-p，则原式$=\frac{{\alpha }^r}{B(a,b)}{\int }_0^1{q}^{b-1}(1-q{)}^{a-1}(\alpha +t-qt{)}^{-r}dq$

$=(\frac{\alpha }{\alpha +t}{)}^r\frac{1}{B(a,b)}{\int }_0^1{q}^{b-1}(1-q{)}^{a-1}(1-\frac{t}{\alpha +t}q{)}^{-r}dq$

而根据超几何函数的积分形式（可以看：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 ）

${}_2{F}_1(a,b;c;z)=\frac{1}{B(b,c-b)}{\int }_0^1{t}^{b-1}(1-t{)}^{c-b-1}(1-zt{)}^{-a}dt,c>b$

原式$=P(在t时刻依然存活|r,\alpha ,a,b)={\left(\frac{\alpha }{\alpha +t}\right)}^r{}_2{F}_1(r,b;a+b;\frac{t}{\alpha +t})$

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

$E(X(t)|\lambda ,p)=\frac{1}{p}-\frac{1}{p}{e}^{-\lambda pt}$

首先让我们求出$E(X(t)|r,\alpha ,p)$。由于$f(\lambda |r,\alpha )=\frac{{\alpha }^r{\lambda }^{r-1}{e}^{-\lambda \alpha }}{\Gamma (r)},\lambda >0$，我们可以求出：

$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\lambda$

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

对于被减数而言,${\int }_0^1\frac{1}{p}\frac{p^{a-1}(1-p{)}^{b-1}}{B(a,b)}dp=\frac{{\int }_0^1{p}^{a-1}(1-p{)}^{b-1}dp}{B(a,b)}=\frac{B(a-1)(b)}{B(a,b)}=\frac{a+b-1}{a-1}$(贝塔函数递推公式)

对于减数而言，${\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^1{p}^{a-2}(1-p{)}^{b-1}(\alpha +pt{)}^{-r}dp$

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

减数=$\begin{array}{rl}& =\frac{{\alpha }^r}{B(a,b)}{\int }_0^1{q}^{b-1}(1-q{)}^{a-2}(\alpha +t-qt{)}^{-r}dq\end{array}$

$={\left(\frac{\alpha }{\alpha +t}\right)}^r\frac{1}{B(a,b)}{\int }_0^1{q}^{b-1}(1-q{)}^{a-2}(1-\frac{t}{\alpha +t}q{)}^{-r}dq$

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

${\left(\frac{\alpha }{\alpha +t}\right)}^r{\frac{B(a-1,b)}{B(a,b)}}_2{F}_1(r,b;a+b-1;\frac{t}{\alpha +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]$

参数求解

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

根据贝塔函数的递推公式，我们有$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-1)=\frac{a}{b+x-1}B(a,b+x)$

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

从而似然函数能够写成$L(r,\alpha ,a,b|X=x,t_x,T)=A_1{A}_2(A_3+{\delta }_{x>0}{A}_4)$的形式

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

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

$A_3={\left(\frac{1}{\alpha +T}\right)}^{r+x}$

$A_4=\left(\frac{a}{b+x-1}\right)(\frac{1}{\alpha +t_x}{)}^{r+x}$

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

$ln(A_1)=ln(\Gamma (r+x))-ln(\Gamma (r))+rln(\alpha )$

$ln(A_2)=ln(\Gamma (a+b))+ln(\Gamma (b+x))-ln(\Gamma (b))-ln(\Gamma (a+b+x))$

$ln(A_3)=-(r+x)ln(\alpha +T)$

$ln(A_4)=ln(a)-ln(b+x-1)-(r+x)ln(\alpha +t_x)$

最终我们有

$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)$

在我给出的那个页面中，作者使用了Excel作为参数求解的工具：

如图所示，B、C、D列3列是我们的原始数据。x代表从第0期到当前期数的总交易次数，$t_x$代表交易期数，注意$t_x=0$时，x也为0，T代表真实时间与设定的第0时刻相差多少。
将4个参数放到B1到B4，然后分别在F到I列求出$A_1$到$A_4$的对数值，然后在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)=x$中的x
而recency要输入输入顾客最近一次购买的时间，即公式中的$t_x$
最后的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,t_x,T)$
针对于以此都没有购买的客户(x=0),我们假定$P(\text{在T时刻依然活跃}T\mid X=0,T,r,\alpha ,a,b)=1$
而对于购买次数为x的客户，在T时刻活跃的概率为：

$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为用户失活的情形，而$(1-p)e^{-\lambda (T-t_x)}$则代表了用户在x次购买之后，在最后一次交易直到T时间没有任何交易，但是依然活跃的概率。

分子分母同时乘以$(1-p{)}^{x-1}{\lambda }^x{e}^{-\lambda t_x}$，

其中分母为：

$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 }$（当x=0时，原式=1，和此处的${\delta }_{x>0}$没有冲突）

而分子=$(1-p{)}^x{\lambda }^x{e}^{-\lambda T}$

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

$P(在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 dp$

而根据贝叶斯公式，我们有：
$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)}$

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

$P(在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 }^x{e}^{-\lambda T}f(\lambda |r,\alpha )f(p|a,b)d\lambda dp$

式中的积分部分可以写为${\int }_0^1(1-p{)}^{x}f(p|a,b)dp{\int }_0^{\infty }{\lambda }^x{e}^{-\lambda T}f(\lambda |r,\alpha )d\lambda$
这2个形式正好对应了我们在求单个客户表达式时使用的“化简式”，于是我们有：

${\int }_0^1{\int }_0^{\infty }(1-p{)}^x{\lambda }^x{e}^{-\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}}$

最终再把$L(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}}$带入，可以得：

$P(在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}}$

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

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

在最初，我们已经推导了$E(Y(t)|\lambda ,p)=\frac{1}{p}-\frac{1}{p}{e}^{-\lambda pt}$

将这个式子乘以$P(在T时刻依然活跃\mid X=x,t_x,T)=\frac{(1-p{)}^x{\lambda }^x{e}^{-\lambda T}}{L(\lambda ,p|X=x,t_x,T)}$

则$E(Y(t)|X=x,t_x,T,\lambda ,p)$

$=\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)}$

$=\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)}$

和上面一样，我们求$E(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 dp$

$=\frac{\mathrm{A}-\mathrm{B}}{L(r,\alpha ,a,b\mid X=x,t_x,T)}$

其中，$A={\int }_0^1{\int }_0^{\infty }{p}^{-1}(1-p{)}^x{\lambda }^x{e}^{-\lambda T}f(\lambda |r,\alpha )f(p|a,b)d\lambda dp$

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

而B$={\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$

$={\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$

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

${\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)$

从而$B=\frac{\Gamma (r+x){\alpha }^r}{\Gamma (r)B(a,b)}{\int }_0^1{p}^{a-2}(1-p{)}^{b+x-1}(\alpha +T+pt{)}^{-(r+x)}dp$

和上面一样，使用$q=1-p$进行还元，得

$\frac{\Gamma (r+x){\alpha }^r}{\Gamma (r)B(a,b)}{\int }_0^1{q}^{b+x-1}(1-q{)}^{a-2}(\alpha +T+t-qt{)}^{-(r+x)}dq$

$=\frac{\Gamma (r+x){\alpha }^r}{\Gamma (r)B(a,b)(\alpha +T+t{)}^{r+x}}{\int }_0^1{q}^{b+x-1}(1-q{)}^{a-2}(1-\frac{t}{\alpha +T+t}q{)}^{-(r+x)}dq$

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

$B=\frac{B(a-1,b+x)}{B(a,b)}{\frac{\Gamma (r+x){\alpha }^r}{\Gamma (r)(\alpha +T+t{)}^{r+x}}}_2{F}_1(r+x,b+x;a+b+x-1;\frac{t}{\alpha +T+t})$

综上，我们可以求得:

$\begin{array}{rl}& 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(r+x,b+x;a+b+x-1;\frac{t}{\alpha +T+t})\right]}{1+{\delta }_{x>0}\frac{a}{b+x-1}{\left(\frac{\alpha +T}{\alpha +t_x}\right)}^{r+x}}\end{array}$

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

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

$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!}]$
而对于用户在T时刻不再活跃的情情形，我们有

$P(Y(t)=y|\lambda ,p)=\{\begin{array}{ll}1& \text{如果}y=0\\ 0& \text{其他情况}\end{array}$

我们分别将$P(在T时刻依然活跃|X=x,t_x,T,\lambda ,p)=\frac{(1-p{)}^x{\lambda }^x{e}^{-\lambda T}}{L(\lambda ,p|X=x,t_x,T)}$

与
$P(在t_x时刻活跃，但是到了第x次交易之后立刻失活)=\frac{p(1-p{)}^{x-1}{\lambda }^x{e}^{-\lambda t_x}}{L(\lambda ,p|X=x,t_x,T)}$

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

$P(Y(t)=y|X=x,t_x,T,\lambda ,p)=\left\{{\delta }_{x>0,y=0}p(1-p{)}^{x-1}{\lambda }^x{e}^{-\lambda t_x}+(1-p{)}^{x+y}\frac{{\lambda }^{x+y}{t}^y{e}^{-\lambda (T+t)}}{y!}+{\delta }_{y>0}p(1-p{)}^{x+y-1}\left[{\lambda }^x{e}^{-\lambda T}-e^{-\lambda (T+t)}{\sum }_{j=0}^{y-1}\frac{{\lambda }^{x+j}{t}^j}{j!}\right]\right\}/L(\lambda ,p|X=x,t_x,T)$

和上面一样，我们需将$\lambda$和p转换为我们求出的那4个参数：

$\begin{array}{rl}P(Y(t)=y|X& =x,t_x,T,r,\alpha ,a,b)\\ & ={\int }_0^1{\int }_0^{\infty }\{P(Y(t)=y|X=x,t_x,T,\lambda ,p)f(\lambda ,p|r,\alpha ,a,b,X=x,t_x,T)\}d\lambda dp\end{array}$

这个$ f(\lambda,p | r,\alpha,a,b,X=x,t_x,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)}$

其中,$A=\frac{B(a+1,b+x-1)}{B(a,b)}\frac{\Gamma (r+x)}{\Gamma (r)}\frac{{\alpha }^r}{(\alpha +t_x{)}^{r+x}}$
,

$B=\frac{B(a,b+x+y)}{B(a,b)}\frac{\Gamma (r+x+y)}{\Gamma (r)y!}\frac{{\alpha }^r{t}^y}{(\alpha +T+t{)}^{r+x+y}}$

$C=\frac{B(a+1,b+x+y-1)}{B(a,b)}\left\{\begin{array}{r}\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 }^r{t}^j}{(\alpha +T+t{)}^{r+x+j}}\end{array}\right\}$

lifetimes中的函数

在lifetimes包中，使用conditional_expected_number_of_purchases_up_to_time函数来预测t期之后的交易数，也就是我们求的$E(Y(t)|X=x,t_x,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,\alpha ,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,y_x,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分析案例。