変分ベイズで RFM 指標から顧客生涯価値を計算してみる

モデルと尤度

RFM 指標から将来の購買回数を予測する Pareto / NBD モデルから派生したモデルに BG / NBD モデルがあります。

Pareto / NBD モデルより計算がかんたんです。

モデル：

顧客の購買はレート $\lambda$ の定常ポアソン過程に従う
顧客の購買のたびに確率 $p$ で生存、 $1-p$ で離脱する
$\lambda$ の事前分布はパラメータ $(\alpha, \beta)$ のガンマ分布とする
$p$ の事前分布はパラメータ $(a,b)$ のベータ分布とする

さらに、 $\lambda$ と $p$ は顧客ひとりひとりで異なるとします。

通常は $\lambda$ と $p$ を周辺化して最尤法でパラメータを推定することが多いようですが、今回は変分ベイズでパラメータを求めます。

変分ベイズを使うメリットは計算の速さです。

記法：

観測期間の終点を $T$ とする
観測期間内の購買回数の合計（フリクエンシー）を $k$ とする
顧客の i 回目の購買日時を $t_i$ とする
$\Delta t_i = t_{i+1}-t_i$ とする

最後の購買時点で顧客が生存している場合を $s=1$ 、離脱している場合を $s=0$ として、考えるべき尤度は、

$L= \prod_{i=1}^{k-1}p\lambda \exp(-\lambda \Delta t_i) \{p\exp(-\lambda (T-t_k))\}^s (1-p)^{1-s}$

です。

対数をとって整理すると、

$\log L = (x-1)\log \lambda +(x-1)\log p- \lambda (t_x-t_1) \\ \quad + s (\log(p)-\lambda(T-t_x)) + (1-s)\log(1-p)$

となります。

変分推論のアルゴリズムの導出

s のサポートが ${0,1}$ であることに注意すると、s の近似事後分布は以下のベルヌーイ分布です。

$q(s) = q^s (1-q)^{1-s}$

q は s にかかっている因子と (1-s) にかかっている因子が足して 1 になるよう正規化したものですから整理すると、以下のようになります。

$\displaystyle q= \frac{\exp(E[\log p])\exp(-E[\lambda] (T-t_k))}{(\exp(E[\log p])\exp(-E[\lambda] (T-t_k)))+\exp(E[\log(1-p)])}\\ \displaystyle= 1/(1+\exp(E[\log(1-p)]/(\exp(E[\log p])\exp(-E[\lambda] (T-t_k))))$

s の期待値は q です。

$\lambda$ の近似事後分布はパラメータ $(k-1, E[s]T+(1-E[s])t_k -t_1)$ のガンマ分布
$\mu$ の近似事後分布はパラメータ $(k-1+E[s]+a, 1-E[s]+b)$ のベータ分布

です。

アルゴリズムをまとめて R で書くとこうなりました。

VIgeom <- function(k,tk,t1,Tim,maxit = 20, alpha=1, beta=1, a=1, b=1){
  N <- length(k)
  p <- runif(N)
  logp <- log(p)
  lambda <- (k+1)/(Tim)
  s <- 1/(1+(1-p)/(p*exp(-lambda*(Tim-tk))))
  for(i in 1:maxit){
    alphahat <- (k-1+alpha)
    betahat <- (s*Tim+(1-s)*tk-t1+beta)
    lambda <- alphahat/betahat
    loglambda <- digamma(alphahat) - log(betahat)
    ahat <- k-1+s+a
    bhat <- 1-s+b
    logp <- digamma(ahat)-digamma(ahat+bhat)
    log1mp <- digamma(bhat)-digamma(ahat+bhat)
    s <- 1/(1+exp(log1mp-logp+lambda*(Tim-tk)))
  }
  list(lambda=lambda,p=ahat/(ahat+bhat),ahat=ahat,bhat=bhat,s=s)
}

生存時間の分布

確率 p で生存する試行を k 回繰り返すので：

$\displaystyle \Pr (\mbox{$T$ 時点で生存}) = \Pr(\mbox{生存時間$>T$}) \\ \displaystyle= \sum_{k=0}^{\infty} p^{k} \frac{(\lambda T)^k \exp(-\lambda T)}{k!}\\ \displaystyle =\exp(-\lambda(1-p)T)$

これは生存関数であるから, 生存時間の密度関数 g(T) は、

$\displaystyle g(T) = \lambda (1-p) \exp(-\lambda(1-p)T)$

です。

有用な統計量

期間 $\tau$ までの期待購買回数は、
$\displaystyle \frac{1}{1-p}(1-\exp(-p\lambda\tau))$

期間 $\tau$ （ $\tau$ 日間）の購買回数を $K(\tau)$ とする。

$\tau$ 時点で生存している場合と $[0,\tau]$ で離脱している場合を考え、
$\displaystyle E[K(\tau)] = \lambda \tau\Pr(T>\tau) + \lambda \int_{0}^{\tau} T g(T) \,dT \\ \displaystyle=\lambda \tau \exp(-\lambda(1-p)\tau) + \int_{0}^{\tau} \lambda T \lambda (1-p) \exp(-\lambda(1-p)T) \,dT\\ \displaystyle=\int^{\tau}_{0}\lambda \exp(-\lambda(1-p)T) \,dT\\ \displaystyle=\frac{1}{1-p}(1-\exp(-\lambda(1-p)\tau))$

（2行目から3行目の変形は部分積分を用いた）

本当はこの統計量の平均を求めたほうが気持ちがいいのですが、めんどうなのでパラメータの平均をプラグインして予測します。

例題

顧客生涯価値(CLV)の計算 with R | かものはしの分析ブログで紹介されている Bruce Hardie: Datasets を用いて将来の購買回数を予測してみます。

タイトルでは顧客生涯価値としましたが、予測するのは単に予測期間の購買回数です。

1997-01-01 から 1997-12-31 までを学習用、1998-01-01から1998-06-30までをテスト用データとします。

library(tidyverse)
library(cowplot)
df <- read_table(file = "~/data/CDNOW_master/CDNOW_master.txt",col_names = FALSE)
df <- mutate(df,X2=as.Date(as.character(X2),"%Y%m%d"))
# range(df$X2)
# "1997-01-01" "1998-06-30"
RFM1 <- dplyr::filter(df,X2<"1998-01-01") %>% 
  arrange(X2) %>% 
  group_by(X1) %>% 
  summarise(first=first(X2),
            last=last(X2),
            frequency = n(),
            monetary=sum(X4)) %>% 
  ungroup() %>% 
  mutate(Tim = as.integer(max(last) - min(first)),
         t1 = as.integer(first - min(first)),
         tx = as.integer(last - min(first)))
RFM2 <- dplyr::filter(df,X2>="1998-01-01") %>% 
  arrange(X2) %>% 
  group_by(X1) %>% 
  summarise(frequency_test = n(),monetary_test=sum(X4))

RFM_join <- left_join(RFM1,RFM2) %>% 
  mutate(frequency_test=if_else(is.na(frequency_test),0L,frequency_test))
outVI <- VIgeom(k = RFM1$frequency,tk = RFM1$tx,t1 = RFM1$t1,Tim = RFM1$Tim,
                alpha = 1,beta=1,a = 1,b=1)

Ti <- as.integer(as.Date("1998-06-30") - as.Date("1998-01-01"))
predfreq <- outVI$s*(1-exp(-outVI$lambda*(1-outVI$p)*Ti))/(1-outVI$p)

compdf <- data.frame(test=RFM_join$frequency_test,predict=predfreq)

ggplot(compdf,aes(x=test,y=predict))+
  geom_hex()+
  geom_abline(slope = 1,intercept = 0,linetype=2)+
  theme_cowplot()

次の図は期待購買回数と実測値の散布図です。

f:id:abrahamcow:20201020040416p:plain

RMSE は1.28でした。

sqrt(mean((RFM_join$frequency_test-predfreq)^2))
# 1.283156

購買回数が0か1以上かという予測が大事だと思ったので、予測値を丸めてアキュラシーを求めてみました。

アキュラシーは約0.8でした。まあまあ妥当な予測と言えそうです。

pred2 <- predfreq>0.5
test2 <- RFM_join$frequency_test>0
tab <- table(pred2,test2)
sum(diag(tab))/sum(tab)
# 0.8071277

ジョンとヨーコのイマジン日記

キョウとアンナのラヴラブダイアイリー改め、ジョンとヨーコのイマジン日記です。

変分ベイズで RFM 指標から顧客生涯価値を計算してみる

モデルと尤度

変分推論のアルゴリズムの導出

生存時間の分布

有用な統計量

例題

おすすめ書籍