読者です 読者をやめる 読者になる 読者になる

廿TT

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

阿部(2004)RF指標から生存確率を求める

R

はじめに

阿部誠(2004)の手法をシミュレーションで再現してみます.
http://merc.e.u-tokyo.ac.jp/mmrc/dp/pdf/MMRC16_2004.pdf

顧客関係管理(CRM)の現場では、RFM 分析という手法が広く使われているそうです. 顧客の購買行動の詳細を蓄積していない企業でも R(リセンシー:直近の購買時点)や F (フリクエンシー:観測期間中の購買回数)だけは把握していることがままあるようです.

データセットはこんな感じです.

       first  recency frequency
1  0.1428815 9.062408        12
2  7.5746912 8.544309         3
3  4.8364402 9.399117        11
4  9.2954182 9.295418         0
5  8.0067610 8.234363         1
6  3.5916135 8.827763        11
7  4.2186687 9.523106        18
8  4.0762378 9.422262        13
9  8.1941630 9.720648         4
10 7.9261058 8.495394         1

first は初回訪問日です.

これは後述するシミュレーションで生成した擬似的なものです.

お客さんの中には途中で離反して店舗に訪れなくなってしまう人もいます.

これを阿部(2004)では死亡と呼び, 顧客の死亡の見極めが大事だとしています.

モデルと尤度関数

阿部(2004)では, RF 指標に対して, 非常にシンプルな仮定を置きます.

仮定 1:お客さんの購買行動はレート λポアソンプロセスに従う.
仮定 2:お客さんの生存時間 τ(購買が維持される期間)はパラメータ μ の指数分布に従う.

この仮定のもとで観測は区間  [0, T_i ] で行われたします.

仮定 1 よりお客さんの購買行動の間隔は指数分布に従うため, 和の分布はガンマ分布になります.

そのため, 顧客 i が時刻  t_i x_i 回目の購買行動を行う密度は,

 \displaystyle \frac{ \lambda ^{x_{i}} t_{i}^{x_i-1}}{\Gamma(x_i)} \exp(-\lambda t_i)

となります.
x_i はフリクエンシー、t_i はリセンシーからわかる)

これより尤度を構成します.

顧客が生存しており, かつx_i 回目の購買が  t_i に起き, かつ区間  [t_i,T_i] に購買が起きない確率は,
\displaystyle e^{-\mu T_i} \times \frac{ \lambda ^{x_{i}} t_{i}^{x_i-1}}{\Gamma(x_i)} \exp(-\lambda t_i) \times  \exp(-\lambda (T_i-t_i)).

(3つめの因子はポアソン分布の確率関数に0を代入した形です. 直近の購買がある日 t だった場合, リセンシーは「時刻 t から観測終了時刻 T までの間に, 購買行動が起こっていない」という情報を持っていることに注意してください.)

顧客が時刻 y_it_i < y_i < T_i )に死亡しており, かつx_i 回目の購買が  t_i に起き, かつ区間  [y_i,T] に購買が起きない確率は,
\displaystyle \int_{t}^{T_i} \mu e^{-\mu y} \times \frac{ \lambda ^{x_{i}} t_{i}^{x_i-1}}{\Gamma(x_i)} \exp(-\lambda t_i) \times \exp(-\lambda (y - t_i)) dy.


生存と死亡という二つの事象は排反なので足し合わせて整理して, 最大化すべき尤度は

\displaystyle L(x_i,t_i,T_i|\lambda,\mu)= \lambda ^{x_i} \left\{ \frac{\lambda}{\lambda+\mu} e^{-(\lambda+\mu)T_i} + \frac{\mu}{\lambda+\mu} e^{-(\lambda+\mu)t_i} \right\}

となります.

シミュレーション

この推定量の精度が気になりますので, R で試してみます.

初回訪問日 o区間 [0, 10] の一様乱数で選びました.

生存時間に対応する値を指数乱数で決め,  \min(o+\tau,10) の間に何回訪問したかを数えるようなプログラムを書きました.

尤度関数の最大化は optim で行いました.

Tim=10
lambda <- 2
mu <- 0.1
simu_rf2 <- function(n=10){
  frequency <- recency <- first <- y <- z <- numeric(n)
  for(i in 1:n){
    o <- runif(1,0,Tim)
    t1 <-cumsum(rexp(1000,lambda))
    tau <- rexp(1,mu) #survival time
    y[i] <- min(o+tau,Tim) 
    tmp <- o + t1<y[i]
    frequency[i] <- sum(tmp)
    first[i] <- ifelse(frequency[i]==0,0,o+t1[1])
    recency[i] <-  ifelse(frequency[i]==0,0,o+t1[frequency[i]])
    z[i] <- as.numeric(o+tau>Tim)
  }
  frequency = frequency-1
  out <- data.frame(first,recency,frequency,y,z)
  out[out$frequency>=0,]
}

objfun <- function(par){ #likelihood
  lambda <- exp(par[1])
  mu <- exp(par[2])
  x <- with(dat,frequency)
  T <- Tim - with(dat,first)
  t <- with(dat,recency) - with(dat,first)
  sum(x*log(lambda) + 
        log( (lambda/(lambda+mu))*exp(-(lambda+mu)*T)+
               (mu/(lambda+mu))*exp(-(lambda+mu)*t) ))
}

res <- vector("list",1000)
for(i in 1:1000){
  dat <-simu_rf2(n=10)
  res[[i]] <- optim(c(1,1),objfun,control=list(fnscale=-1))
}
library(rlist)
table(unlist(list.map(res,convergence)))  #check convergence
lambdahat <- exp(unlist(list.map(res,par[1])))
muhat <- exp(unlist(list.map(res,par[2])))
hist(lambdahat,breaks="scott")
abline(v=mean(lambdahat),col="royalblue",lwd=2)
abline(v=lambda,col="tomato",lwd=2) #true
hist(muhat,breaks="scott")
abline(v=mean(muhat),col="royalblue",lwd=2)
abline(v=mu,col="tomato",lwd=2) #true

シミュレーションの試行回数を1000回として, サンプルサイズ10のとき, λ の推定値は以下のヒストグラムのように分布しています.

f:id:abrahamcow:20160405043828p:plain

赤い線がシミュレーションで仮定した真値, 青い線が推定値の平均です.

同様に μ の推定値のヒストグラムは以下です.

f:id:abrahamcow:20160405044108p:plain

若干正のバイアスがあるようですが, 概ね真値に近い値が求まっています.

阿部(2004)では, パラメータ λμ は顧客ごとに異なるとして, 顧客ひとりひとりにたいして, 推定を行っています.

res2 <- vector("list",1000)
for(i in 1:1000){
  dat <-simu_rf2(n=1)
  res2[[i]] <- optim(c(1,1),objfun,control=list(fnscale=-1))
}
table(unlist(list.map(res2,convergence)))  #check convergence
lambdahat <- exp(unlist(list.map(res2,par[1])))
muhat <- exp(unlist(list.map(res2,par[2])))
hist(lambdahat,breaks="scott")
abline(v=mean(lambdahat),col="royalblue",lwd=2)
abline(v=lambda,col="tomato",lwd=2) #true
hist(muhat,breaks="scott")
abline(v=mean(muhat),col="royalblue",lwd=2)
abline(v=mu,col="tomato",lwd=2) #true

サンプルサイズ1のとき, λμ の推定値は以下のヒストグラムのように分布しています.

f:id:abrahamcow:20160405044557p:plain

f:id:abrahamcow:20160405044638p:plain

μ の推定値はかなりバラつきも, バイアスも大きくなってしまうようです.

阿部はベイズ推定を用いることで, この課題を解決できるとしていますが, 今のところベイズ推定はぼくの能力不足で再現できませんでした.

どなたかやってみてください.