1. BG/NBD概率模型
BG/NBD 模型又称为贝塔几何/负二项模型。他是基于 Pareto/NBD 模型假设设计的概率预测模型。BG/NBD 模型是用于描述非契约客户关系情境下重复购买行为。即用户可以随时购买产品，无时间约束。该模型可利用用户历史交易数据 (RFM) 来预测未来每个用户的交易次数和流失率， BG/NBD 模型主要有以下五个假设：
假设1：用户在活跃状态下，一个用户在时间段 $t$ 内完成的交易数量服从均值为 $λt$ 的泊松分布。泊松分布的 PDF 函数如下，$j$ 为交易次数。比如一个顾客在 $λ=5$ 时的概率分布图如下所示：
$$f(t_j|t_{j-1};\lambda)=\lambda e^{-\lambda(t_j-t_{j-1})}，t_j>t_{j-1} \geq 0$$
import matplotlib.pyplot as plt
from scipy.stats import poisson
probability_arr = []
distribution = poisson(5)
for transactions in range(0,20):
probability_arr.append(distribution.pmf(transactions))
plt.figure(figsize=(8,5))
plt.ylabel('Probability')
plt.xlabel('Number of Transactions')
plt.xticks(range(0, 20))
plt.title('Probability Distribution Curve')
plt.plot(probability_arr, color='black', linewidth=0.7, zorder=1)
plt.scatter(range(0, 20), probability_arr, color='purple', edgecolor='black', linewidth=0.7, zorder=2)
plt.show()
图 1 顾客在λ=5时的概率分布图
假设2：用户的交易率 $λ$ 服从形状参数为 $r$，逆尺度参数为 $α$ 的 gamma 分布，PDF 函数如下所示。比如 100 个用户的交易率 $λ$ 服从 $r=9.0，α=0.5$ 的gamma分布时，这 100 个用户的泊松分布图如下所示：
$$f(\lambda|r,\alpha)=\frac{\alpha^r\lambda^{r-1}e^{-\lambda\alpha}}{\tau(r)}, \lambda>0$$
import numpy as np
plt.figure(figsize=(8,5))
for customer in range(0, 100):
distribution = poisson(np.random.gamma(shape=9, scale=0.5))
probability_arr = []
for transactions in range(0,20):
probability_arr.append(distribution.pmf(transactions))
plt.plot(probability_arr, color='black', linewidth=0.7, zorder=1)
plt.ylabel('Probability')
plt.xlabel('Number of Transactions')
plt.xticks(range(0, 20))
plt.title('Probability Distribution Curve 100 Customers')
plt.show()
图 2 交易率 λ 服从 r=9.0，α=0.5 的 gamma 分布时的泊松分布图
假设3：每个用户在交易 $j$ 完成后流失的概率服从参数为 $p$(流失率)的几何分布，PDF 函数如下所示。
$$P(第 j 次交易后的流失率)=p(1-p)^{j-1},j=1,2,3,...$$
假设4：用户的流失率 $p$ 服从形状参数为 $a,b$ 的 beta 分布，PDF 函数如下所示。比如 100 个用户的流失率 $p$ 服从 $a=1.0，b=2.5$ 的 beta 分布时，这 100 个用户的几何分布图如下所示：
$$f(p|a,b)=\frac{p^{a-1}(1-p)^{b-1}}{B(a,b)}, 1 \geq p \geq 0$$
import numpy as np
plt.figure(figsize=(8,5))
for customer in range(0, 100):
distribution = poisson(np.random.gamma(shape=9, scale=0.5))
probability_arr = []
beta = np.random.beta(a=1.0, b=2.5)
cumulative_beta = 0
for transactions in range(0,20):
proba = distribution.pmf(transactions)
cumulative_beta = beta + cumulative_beta - (beta * cumulative_beta)
inactive_probability = 1 - cumulative_beta
proba *= inactive_probability
probability_arr.append(proba)
probability_arr = np.array(probability_arr)
probability_arr /= probability_arr.sum()
plt.plot(probability_arr, color='black', linewidth=0.7, zorder=1)
plt.ylabel('Probability')
plt.xlabel('Number of Transactions')
plt.xticks(range(0, 20))
plt.title('Probability Distribution Curve 100 Customers with drop-off probability after each transaction')
plt.show()
图3 用户交易流失率
假设5：每个用户的交易率 $λ$ 和流失率 $p$ 互相独立。
2. 用户生命周期预测，Python实现
图 5 用户交易生命周期
图 5 数据展示
模型原理：首先 gamma 分布和 beta 分布分别为参数交易率 $λ$ 和流失率 $p$ 的先验分布，而泊松分布和几何分布是样本的分布函数，即似然函数。接下来建立交易率 $λ$ 和流失率 $p$ 的联立似然函数，使用 Nelder-Mead 的单纯形算法求解 gamma 分布和 beta 分布中的参数 (r, α, a, b)，这是一种启发式的，非梯度搜索方法来最小化负对数似然代价函数。
$$\ln[L(r,\alpha,a,b|X=x,t_x,T)=\ln(A_1)\ln(A_2)\ln(e^{ln(A_3)}+\delta_{x>0}e^{\ln(A_4)})],\,\,where$$
$$\ln(A_1)=\ln[\tau(r+x)]-\ln[\tau(r)]+r\ln(\alpha)$$
$$\ln(A_2)=\ln[\tau(a+b)]+\ln[\tau(b+x)]-\ln[(\tau(b)]-\ln[\tau(a+b+x)]$$
$$\ln(A_3)=-(r+x)\ln(\alpha+T)$$
$$\ln(A_4)=\begin{cases} \ln(a)-\ln(b+x-1)-(r+x)\ln(\alpha+t_x), if\,\,x>0\\ 0,\,\,otherwise\end{cases} $$
用代码实现上面的似然函数如下：
from scipy.special import gammaln
def negative_log_likelihood(params, x, t_x, T):
if np.any(np.asarray(params) <= 0):
return np.inf
r, alpha, a, b = params
ln_A_1 = gammaln(r + x) - gammaln(r) + r * np.log(alpha)
ln_A_2 = (gammaln(a + b) + gammaln(b + x) - gammaln(b) -
gammaln(a + b + x))
ln_A_3 = -(r + x) * np.log(alpha + T)
ln_A_4 = x.copy()
ln_A_4[ln_A_4 > 0] = (
np.log(a) -
np.log(b + ln_A_4[ln_A_4 > 0] - 1) -
(r + ln_A_4[ln_A_4 > 0]) * np.log(alpha + t_x)
)
delta = np.where(x>0, 1, 0)
log_likelihood = ln_A_1 + ln_A_2 + np.log(np.exp(ln_A_3) + delta * np.exp(ln_A_4))
return -log_likelihood.sum()
参数求解代码如下：
from scipy.optimize import minimize
scale = 1 / df['T'].max()
scaled_recency = df['t_x'] * scale
scaled_T = df['T'] * scale
def _func_caller(params, func_args, function):
return function(params, *func_args)
current_init_params = np.array([1.0, 1.0, 1.0, 1.0])
output = minimize(
_func_caller,
method="Nelder-Mead",
tol=0.0001,
x0=current_init_params,
args=([df['x'], scaled_recency, scaled_T], negative_log_likelihood),
options={'maxiter': 2000}
)
r = output.x[0]
alpha = output.x[1]
a = output.x[2]
b = output.x[3]
alpha /= scale
print("r = {}".format(r))
print("alpha = {}".format(alpha))
print("a = {}".format(a))
print("b = {}".format(b))
计算得到r=0.243，α=4.414，a=0.793，b=2.426.
用户总交易次数预测
接下来通过使用上面的四个参数建立预测模型，即求解交易次数的期望 E(x)。$_2F_1$ 为高斯超几何函数。
$$E(X(t)|r,\alpha,a,b)=\frac{a+b-1}{a-1}[1-(\frac{\alpha}{\alpha+t})^r_2F_1(r,b;a+b-1;\frac{t}{\alpha+t})]$$
$_2F_1$为高斯超参函数，函数形式为：
$$_2F_1(a,b;c;z)=\sum_{j=0}^{\infty}\mu_j ,\,\,where \,\, \mu_j=\frac{(a)_j}{(b)_j}\frac{z^j}{j!}$$
上面的公式可以展开为：
$$\frac{\mu_j}{\mu_{j-1}}=\frac{(a+j-1)(b+j-1)}{(c+j-1)j}z$$
from scipy.special import hyp2f1
def expected_sales_to_time_t(t):
hyp2f1_a = r
hyp2f1_b = b
hyp2f1_c = a + b - 1
hyp2f1_z = t / (alpha + t)
hyp_term = hyp2f1(hyp2f1_a, hyp2f1_b, hyp2f1_c, hyp2f1_z)
return ((a + b - 1) / (a - 1)) * (1-(((alpha / (alpha+t)) ** r) * hyp_term))
expected_sales_to_time_t(78)
现在通过以上代码预测未来 78 周单个用户的交易次数为 1.858 次，但计算 E(x) 为总的用户购买总次数，这里不能简单的通过单个用户交易次数乘以总用户数得到，因为每个用户的第一次交易时间点不同。这里统计有 ns 个用户在第 s 天进行了第一次交易。那么这 ns 个用户在未来某段时间内的交易次数是相同的。
$$t=\sum^S_{s=0}\delta_{(t>s)}n_sE[X(t-s)],where\,\, \delta_{(t>s)} \begin{cases} 1,\,\,if\,\,t>s \\ 0,\,\,otherwise \end{cases}$$
# Period of consideration is 39 weeks. # T indicates the length of time since first purchase n_s = (39 - df['T']).value_counts().sort_index() n_s.head()
数据得到 18 个用户在第 0.14 周进行了第一次交易，22 个用户在第 0.285 周进行了第一次交易，17 个用户在第 0.428 周进行了第一次交易。那么接下来只需要计算每个交易时间下的未来交易量再乘以相同交易时间下的用户数，最后求和即可得到总的交易次数。
## 1/7.0 is a day.
forecast_range = np.arange(0, 78, 1/7.0)
def cumulative_repeat_transactions_to_t(t):
expected_transactions_per_customer = (t - n_s.index).map(lambda x: expected_sales_to_time_t(x) if x > 0 else 0)
expected_transactions_all_customers = (expected_transactions_per_customer * n_s).values
return expected_transactions_all_customers.sum()
cum_rpt_sales = pd.Series(map(cumulative_repeat_transactions_to_t, forecast_range), index=forecast_range)
cum_rpt_sales.tail(10)
最后算出接下来 78 周用户的交易总次数为 4156 次。通过了解到未来用户的交易次数，可以计算未来的收入，从而在现阶段，运营人员可以更好地计算推广预算制定相应的运营推广方案。
每个用户交易次数的条件预测
为预测每个用户在未来一段时间内的交易次数，这里推导出条件期望，根据用户历史的交易次数和交易时间数据，并根据上面得到的分布函数参数值，条件期望的最终计算公式如下所示：
$$E(Y(t)|X=x,t_x,T,r,\alpha,a,b)=\frac{a+b+x-1}{a-1}*\frac{[1-(\frac{\alpha+T}{\alpha+T+t})^{r+x}_2F_1(r+x,b+x;a+b+x-1;\frac{t}{\alpha+T+t})]}{1+\delta_{x>0}\frac{a}{b+x-1}(\frac{\alpha+T}{\alpha+t_x})^{r+x}}$$
def calculate_conditional_expectation(t, x, t_x, T):
first_term = (a + b + x - 1) / (a-1)
hyp2f1_a = r + x
hyp2f1_b = b + x
hyp2f1_c = a + b + x - 1
hyp2f1_z = t / (alpha + T + t)
hyp_term = hyp2f1(hyp2f1_a, hyp2f1_b, hyp2f1_c, hyp2f1_z)
second_term = (1 - ((alpha + T) / (alpha + T + t)) ** (r + x) * hyp_term)
delta = 1 if x > 0 else 0
denominator = 1 + delta * (a / (b + x - 1)) * ((alpha + T) / (alpha + t_x)) ** (r+x)
return first_term * second_term / denominator
calculate_conditional_expectation(39,2,30.43,38.86)
通过上面的函数代入某用户过去 38.86 周内 (T=38.86) 交易两次 (x=2)，第二次交易时间为 30.43 周(tx=30.43) 的条件，计算得到该用户在未来 39 周将交易 1.226 次。通过预测出每个用户未来的交易次数，可以更针对性地细分用户人群，找准目标价值人群从而制定细分运营方案，比如未来一年 52 周用户预测出将交易 1 次，那么该用户有流失的风险，那么在现阶段实施促销方案（如发放促销卡），提高用户的交易频次将减小用户流失的风险。
4. 学习资料
- 论文推荐：代码中都直接计算最终的预测公式，公式的推导过程可详细查看该论文。http://brucehardie.com/papers/018/fader_et_al_mksc_05.pdf
- Python语言lifetimes库：https://pypi.org/project/Lifetimes/
- 参考博客：https://benalexkeen.com/bg-nbd-model-for-customer-base-analysis-in-python/
评论