廿TT

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

ポアソン過程・フラットモデル(BTYDモデル番外編)

今日の川柳

モデルと尤度

RFM 指標から将来の購買回数を予測するモデルとして、次の2つを紹介しました。

ポイントは観測されない顧客の生存と離脱を潜在変数としてモデルに含めるところです。

このモデルをもうちょっとシンプルにできないかなと思って以下のモデルを試してみました。

モデル:

  1. 顧客の購買はレート  \lambda の定常ポアソン過程に従う
  2. 顧客の生存時間は区間[0,\infty) の一様分布に従う
  3. \lambda\mu の事前分布はパラメータ (\alpha,\beta) のガンマ分布とする

 \lambda は顧客ひとりひとりで異なるとします。

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

記法:

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

考えるべき尤度は、

L= \prod_{i=1}^{x-1}\lambda \exp(-\lambda \Delta t_i) \{\exp(-\lambda (T-t_x))\}^(I_{y>T}) \{\exp(-\lambda (y-t_x))\}^{I_{t_x < y \le T}}

です。

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

 \log L = (x-1)\log \lambda - \lambda (t_x-t_1) + I_{y>T} (-\lambda(T-t_x)) + I_{t_x < y \le T} ( -\lambda (y-t_x))

となります。

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

まず  \lambda を所与として、 y の近似事後分布の正規化定数を求めます。

式の展開では省略しますが、以下で単に \lambda となっているところは、正確には  E[\lambda] であることに注意してください。

y に依存しない項を省略して書くと

 q(y) \propto \exp( -\lambda y)^{I_{t_x < y \le T}}

です。これを積分して

 \int^{T}_{t_x} \exp( -(\lambda)y) dy=[ -\frac{1}{\lambda}\exp( -\lambda y)]^{T}_{t_x} \\
=\frac{1}{\lambda} \{ \exp( -\lambda t_x) - \exp( -\lambda T)

これを C とおくと

 E[y] = \frac{1}{C} \int^T_{t_x} y \exp(-(\lambda)y) dy\\
=  \frac{1}{C(\lambda)^2} [-\exp(-(\lambda))\{(\lambda)y+1\} ]^{T}_{t_x}\\
= \frac{1}{\lambda}\frac{\exp(-(\lambda)t_x)\{(\lambda)t_x+1\}-\exp(-(\lambda)T)(\{\lambda\}T+1)}{
      \exp(-\lambda t_x)-\exp(- \lambda T)}

です。

次に、s= I_{y>T} とおいて s の近似事後分布をもとめます。

y と s が独立として近似することも可能ですが、さすがに無理があると思うので、y を周辺化した場合の s の近似事後分布をもとめることにします。

 \int^{T}_{t_x} \exp (- E[\lambda](y- t_x)) =  \frac{1}{\lambda}\{1-\exp(-\lambda(T-t_x))

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

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

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

 p= 1/(1+ (\exp(\{E[\lambda])(T-t_x)\}-1)/E[\lambda])

s の期待値は p です。

さらに、  \lambda の近似事後分布はパラメータ (x-1+\alpha, E[sT+(1-E[s])E[y] -t_1)+\beta] のガンマ分布です。

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

doVI <-function(x,tx,t1,Tim, a_lam =1, b_lam = 1){
  lambda <- (x+1)/Tim
  for (i in 1:20) {
    u <- (1/lambda)*(exp(-lambda*tx)*(lambda*tx+1)-exp(-lambda*Tim)*(lambda*Tim+1))/
      (exp(-lambda*tx)-exp(-lambda*Tim))
    s <- 1/(1+(exp(lambda*(Tim-tx))-1)/lambda)
    lambda<- (x-1+a_lam)/(s*Tim+(1-s)*u-t1+b_lam)
  }
  list(lambda=lambda,s=s,u=u)
}

例題

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

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

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

次の図はROCカーブです。

f:id:abrahamcow:20190813223631p:plain

outVI <- doVI(RFM1$frequency,RFM1$tx,RFM1$t1,RFM1$Tim,a_lam = 2,b_lam = 1)
uid <-unique(RFM2$X1)

Ti <- as.integer(as.Date(max(df$X2)) - as.Date("1998-01-01"))

pred <- prediction(outVI$s*(1-exp(-outVI$lambda*Ti)),(RFM1$X1 %in% uid))
perf <- performance(pred, "tpr", "fpr")
plot(perf)
perf <- performance(pred, "auc")
perf@y.values

AUCは約0.79でした。

まあまあ妥当な予測と言えそうです。

WikipediaBuy Till you Die - Wikipedia)によると他にも顧客の購買をベルヌーイプロセスで表すモデルなどがあるようです。

おすすめ書籍

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