廿TT

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

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

今日の川柳

モデルと尤度

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

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

モデル:

  1. 顧客の購買はレート  \lambda の定常ポアソン過程に従う
  2. 顧客の生存時間はレート  \mu の指数分布に従う
  3. \lambda\mu の事前分布はそれぞれパラメータ  (\alpha,\beta) (a,b) のガンマ分布とする

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

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

記法:

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

ポアソン過程では、顧客の購買間隔は指数分布に従うので、購買間隔に関する尤度は
 \prod_{i=1}^{x-1}\lambda \exp(-\lambda \Delta t_i)
となるように思えますが、観測期間内の最後の購買日(リセンシー)から観測期間の終点までの期間は、「この期間は購買がなかった」という情報を持っているので、これも尤度にかけなければなりません。

I_A を指標関数として、顧客が T 時点で生存している場合と、期間 (t_x,T] で離脱している場合を分けて考えると、購買間隔に関する尤度は

 \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_{y>T}}

となります。新しい因子は指数分布の生存関数です。

y が出てきてしまったので、さらに生存時間 y についての尤度をかけ合わせる必要があるため、最終的に考えるべき尤度は、

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

です。

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

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

となります。

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

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

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

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

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

です。これを積分して

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

これを C とおくと

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

です。

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

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

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

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

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

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

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

s の期待値は p です。

ここまでくればあとはかんたんです。

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

です。

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

doVI <-function(x,tx,t1,Tim,a_lam = 1, b_lam = 1, a_mu = 1, b_mu = 1){
  lambda <- ((x+1)/Tim)
  mu <- 1/(tx+1)
  elmu<-mu
  n <- length(x)
  for (i in 1:20) {
    u <- ifelse(Tim==tx,0,(1/(lambda+mu))*(exp(-(lambda+mu)*(tx))*((lambda+mu)*(tx)+1)-exp(-(lambda+mu)*(Tim))*((lambda+mu)*(Tim)+1))/
      (exp(-(lambda+mu)*(tx))-exp(-(lambda+mu)*(Tim))))
    s <- 1/(1+elmu*(exp((lambda+mu)*(Tim-tx))-1)/(lambda+mu))
    lambda <- ((x-1)+a_lam)/((s*(Tim)+(1-s)*u-t1)+b_lam)
    mu <- ((1-s)+a_mu)/((s*Tim+(1-s)*u)+b_mu)
    elmu <- exp(digamma((1-s)+a_mu)-log((s*Tim+(1-s)*u)+b_mu))
  }
  list(lambda=lambda,mu=mu,s=s,u=u)
}

有用な統計量

期間 \tau までの期待購買回数は、
 (\lambda/\mu)(1-\exp(-\mu \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:20190812232645p: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.289151

RMSE は約1.23 でした。

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

おすすめ書籍

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