干过一些流水预测的活儿
0x00 Introduction
前一阵在某个手游厂做了一些关于流水预测的工作，这里记录一下。
流水预测，这个课题字面意思就是为了预测后续一段时间内公司的流水情况，比如大到一个商场，小到一个小区门口的小卖部，都会想知道自己未来比如1个月的流水如何。
0x01 Keywords
CLV(customer lifetime value) ; LTV(lifetime value) ; BTYD(Buy Till You Die)
0x02 Problem
把流水分成两部分，一部分来自已有用户的贡献，另一部分是目前未知的新用户带来的流水。
这里讨论的只是对于老玩家流水的预测。
在对老玩家LTV预测的问题上，一个分类维度是是否有合同：Contractual/Non-Contractual。有合同指的是我能知道我当前有那些用户，他们分别会在什么时候流失，即合同终止时间，这个适用于比如会员制的超市。实际情况下，对于互联网时代，更多的是Non-Contractual情况，并不能确定的知道用户/玩家什么时候就流失了。
这里我所遇到的也是 Non-contractual Problem。
0x03 Approaches
1. Naive Method
在手游这个domain，比如要预测未来一天的流水，通常人们采取的做法是：根据一些历史经验，估计出那天的DAU，然后乘上付费玩家占比，再乘上他们的平均付费值
2. BTYD
1987年，MIT的David C. Schmittlein发了一篇文章【1】其中提出了一个生成模型，用来预测未来一段时间内老用户的购买次数。模型称作PNBD (Pareto/NBD)。
这篇文章的Intuition是根据每个用户历史购买信息，计算出这个用户在T时刻的存活概率 P(Alive | Information) * 这个用户在(0,T] 的总计购买次数期望 = 这个用户从出生(0时刻)到T的累计购买次数预测。
这里David基于5条假设来给出他的模型：
对于每个用户个体层面:
- 每个用户当还存活时的购买过程服从Poisson Process，参数为$\lambda$ ($1/\lambda$代表两次购买的平均间隔)
- 每个用户随着时间推移，还存活的概率服从Exponential distribution，参数为$\mu$ (代表死亡率)
在整体层面:
- 所有用户的购买过程参数$\lambda$整体服从Gamma distribution: $\Gamma(r,\alpha)$
- 所有用户的死亡率参数$\mu$整体服从另一个Gamma distribution: $\Gamma(s,\beta)$
- $\lambda$与$\mu$独立
在David的推导中，他所用的输入是用户的过去购买时间次数特征。
上图表示了某一个用户在(0,T]时间内一共有x次购买，第一次发生在t1时刻。其中我们把最后一次购买时间tx和第一次购买时间t1之间的时间间隔极为recency（这里时间用的是天为单位，因情况不同而异）。从t1时刻到需要观测的节点T之间的时间间隔，计作t。
这里T为观测结束的时间，“0”表示这个人的“诞生”时刻。在传统的商业模式下（David老爷子的paper可是在1987年的环境下）用户第一次来买之前公司不知道他的存在，所以第一次购买时间作为0时刻，即看作t0。
我们这里对于“诞生”的定义，结合手游场景，我们尝试了用户注册以及第一次购买分别作为“0”时刻，差别不大，最后和原文一样选取第一次购买作为0时刻（t1因此变成第二次购买）。
PNBD就是根据recency, t, T这三个指标作为输入Information进行拟合（训练），在模型拟合完毕之后，可以对于任意一个用户后续任意一天进行购买次数的预测。
>>>>>
PNBD细节原理+工程实践（这里的公式只是为了备忘，可读性比较差，建议跳过或者根据原文推导）：
在David老爷子的假设下，
Naive的想法：对于个体层面，既然每个用户的购买服从Poisson Process，存活率服从Exponential Distribution，那我们用每个人的Information（recency，t，T），对每个用户去拟合这俩分布的参数$\lambda$和$\mu$，然后再把这个分布用来预测。但是这样会有问题，很多用户的购买记录可能非常稀疏（的确很多人只买过一次），直接拟合这俩参数很不靠谱。
原文的做法：既然还有整体层面的3条假设，那我们就不直接拟合$\lambda$和$\mu$。而是转向拿所有人的数据去拟合两个Gamma Distribution的系数。
公式推导上【1】的方式比较复杂，分情况进行了讨论，把用户的实际流失时间$\tau$在T的前面和T之后进行的推导。【2】中在推导上进行了改进，用更简单的表达式整理了结果。并且【8】的计算代码也是根据【2】进行实现的
在省略了亿步推导过后，可得：也是【2】中的式(41)，即：
上式表示在[T,T+t) 的时间段内，任意一个用户的预计购买次数。
其中:
$L(r,\alpha,s,\beta | x,t_x,T) = \frac{\Gamma(r+x)\alpha^r\beta^s}{\Gamma(r)}\{\frac{1}{(\alpha+T)^{r+x}(\beta+T)^s}+(\frac{s}{r+s+x})A_0\}$
当$\alpha \geq \beta : $
$\begin{align} A_0 = &\; \frac{{}_2F_1(r+s+x,s+1;r+s+x+1;\frac{\alpha-\beta}{\alpha+t_x})}{(\alpha+t_x)^{r+s+x}} \\ &- \frac{{}_2F_1(r+s+x,s+1;r+s+x+1;\frac{\alpha-\beta}{\alpha+T})}{(\alpha+T)^{r+s+x}}\end{align}$
当$\alpha \leq \beta : $
$\begin{align} A_0 = &\; \frac{{}_2F_1(r+s+x,r+x;r+s+x+1;\frac{\beta-\alpha}{\beta+t_x})}{(\beta+t_x)^{r+s+x}} \\ &- \frac{{}_2F_1(r+s+x,r+x;r+s+x+1;\frac{\beta-\alpha}{\beta+T})}{(\beta+T)^{r+s+x}}\end{align}$
上面是在参数$r,\alpha,s,\beta $拟合完毕以后用于预测的式子，至于这四个参数拟合，【1】中是靠最大似然估计得到，即目标函数就是Maximize下式：
$LL(r,\alpha,s,\beta) = \sum\limits_{i=1}^{N} \ln[L(r,\alpha,s,\beta | x_i,t_{x_i},T_i)]$
参数估计 -> Optimization Problem
数学上这个问题于是变成了一个优化问题。
在【8】中，求解优化问题的实现用的是scipy的optimize包来实现的，具体scipy的文档可以看【9】。优化问题是个大坑，这里浅尝辄止，记录一下所了解的。
【8】的代码中默认调用的方法是Scipy中的Nelder-Mead方法，这是一种单纯形法。优化问题的求解在一个维度上可以分成两类：Gradient-free / Gradient-involved 即是否需要用到导数信息，单纯形法就属于gradient-free的方法。（大家在DL中耳熟能详的梯度下降，牛顿法都是Gradient-involved求解方法）
两者使用场景不同，比如当目标函数不可导的时候，不需要梯度就是一个很大的优势。当然代价就是速度慢。这里的数据量在100w左右，一个用户的information就是1条，直接丢进【8】的实现中（Nelder-Mead），参数拟合会非常耗时，于是我想着尝试其他求解方式。
除了Nelder-Mead之外Scipy包里的其他的方式基本都是Gradient-Involved。
那问题来了，目标函数的导数怎么计算呢？作为一个数学渣我尝试写出这个目标函数的导数。然而在目标函数中见到${}_2F_1()$这个高斯超几何函数的时候，我心灰意冷，经查，这函数导数没有解析解。
这就意味着不能用gradient-involved的求解方法了吗？也不是，可以通过原始的有限差分法来求解需要的点的导数值: $ y' = \frac{y(x+\epsilon) - y(x)}{\epsilon}$
我通过采样把数据量固定在1w，我运行了5次Nelder-Mead算法，所需要的时间在5～13秒（时间会因为初始值选取而波动）
之后尝试了【9】中scipy的其他所有的方法以后(Powell, Cobyla, SLSQP, CG, BFGS, TNC, L-BFGS-B)，结果SLSQP表现最好，5次均在1～2秒完成。放回到100w+的数据量上也是10分钟内能拟合完成。
当然，SLSQP的代价又是什么呢？在某些随机初始值的选取下会不收敛。
不收敛的问题在5次1w数据量的测试中并没有出现，但是随着数据量变成100w，这个问题显现了，而且出现频率有20%左右，是一个不能忽视的问题。
不收敛 -> 多起始点求解
scipy的optimize包在用SLSQP方式求解的时候是单线程的，cpu0有难，7核围观。这里我对【8】进行了多线程的改造，在参数拟合的时候选取7个不同随机起始点进行迭代，当有一个线程（收敛）找到了最优解以后，就停止其他线程的迭代。这样就把不收敛的概率从0.2降低为 ${0.2}^7$, 即有0.9999872的概率可以成功收敛。
<<<<<
在有办法得到每个用户在后续一段时间内的购买次数以后，我们需要引入金额来实现流水预测。
一个非常拍脑袋的做法就是把每个用户的历史购买平均单价作为这个人的均单价，这好吗，这不好。
比起用样本均值，刻画出历史购买单价的分布，再根据分布得到每个用户的均单价期望更准确。
Peter 在2004年提出了一种计算每个用户均单价的方式【4】，同样他对于回归模型嗤之以鼻，也提出了一个概率模型来刻画用户的均单价。模型叫做GG (Gamma-Gamma)
类似PNBD，GG模型也建立在作者Peter提出的5个假设之上：
对于每个个体：
- 每个客户每笔订单的价格围绕着他的均单价随机分布
- 均单价因客户而异，每个人的均单价在时序上不变
- 每个用户的均单价与其订单数量/购买过程无关
- 每个用户的历史购买单价$z_i$ ~ $\Gamma(p,v)$
- $v$ ~$ \Gamma(q,\gamma) $
>>>>>
同样这里充满公式，可读性很弱，在知道整体是干嘛的以后，跳过即可，想看推导的可以在【4】【5】得到满足。
可以看出这里第4，5条的设定也避免了对于每个个体用户参数的估计，避免了稀疏个体的问题。
在省略了亿步推导之后，得到下式，即【5】中的(5):
$E(Z|p,q,\gamma;\tilde{z}) = (\frac{q-1}{px+q-1})\frac{p\gamma}{q-1} + (\frac{px}{px+q-1})\tilde{z}$
上式代表着任意一个用户的均单价期望。（*即模型predict）
其中$p,q,\gamma$这三个参数估计也同样通过最大似然估计得到
即Maximize下式: (*即模型拟合)
$LL(p,q,\gamma|data) = \sum\limits_{i=1}^N \ln[f(\tilde{z}_i | p,q,\gamma;x_i)]$
其中:
$f(\tilde{z} | p,q,\gamma;x) = \frac{\Gamma(px+q)}{\Gamma(px)\Gamma(q)} \frac{\tilde{z}^{px-1}x^{px}\gamma^q}{(\gamma+x\tilde{z})^{px+q}}$
<<<<<
好了，我们费了很大功夫读 paper，看源码，魔改，兼顾了模型拟合和执行效率，终于把BTYD这一套在流水预测领域里的硬骨头啃了下来，在我们自己的数据集上成功复现。
这好吗，这很好 感谢兄贵分享
ReplyDelete