廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

BG / NBD モデルの変分推論(RFM指標から将来の購買回数を予測する)

今日の川柳

モデルと尤度

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

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

モデル:

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

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

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

記法:

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

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

L= \prod_{i=1}^{x-1}p\lambda \exp(-\lambda \Delta t_i) \{p\exp(-\lambda (T-t_x))\}^s (1-p)^{1-s}

です。

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

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

となります。

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

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

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

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

 p= 1/(1+\exp(E[\log(1-p)]/(\exp(E[\log p])\exp(-\lambda (T-t_k))))

s の期待値は p です。

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

です。

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

doVI <- 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 <- rexp(N)
  theta <- rexp(1)
  logtheta <- log(theta)
  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)/(exp(logp)*exp(-lambda*(Tim-tk))))
  }
  list(lambda=lambda,p=ahat/(ahat+bhat),alpha=alpha,s=s)
}

有用な統計量

期間 \tau までの期待購買回数は、
 \frac{1}{1-p}(1-\exp(-p\lambda\tau))
(計算略: \tau 時点で生存している場合と離脱している場合のふたつを足し合わせ、部分積分を用いれば求まる)

例題

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

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

library(tidyverse)
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()) %>% 
  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())

RFM_join <- left_join(RFM1,RFM2) %>% 
  mutate(frequency_test=if_else(is.na(frequency_test),0L,frequency_test))

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

f:id:abrahamcow:20190813074531p:plain

outVI <- doVI(RFM1$frequency,RFM1$tx,RFM1$t1,RFM1$Tim)
predfreq <- outVI$s*(outVI$lambda/outVI$mu)*(1-exp(-outVI$mu*Ti))

plot(predfreq,RFM_join$frequency_test)
abline(0,1)

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

RMSE は約1.28 でした。

計算も楽だし、予測精度もほとんど変わらないようなので、Pareto / NBD モデルを使う場面で BG / NBD モデルを使っても問題なさそうです。

Pareto / NBD モデルの変分推論(RFM指標から将来の購買回数を予測する) - 廿TT

おすすめ書籍

変分ベイズの解説としては以下の本がわかりやすかったです。