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

廿TT

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

RF指標から平均購買頻度を求める

はじめに

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

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

    recency frequency
1  9.089831         3
2  9.739291         6
3  8.447093         4
4  9.501840         3
5  8.000726         3
6  3.113179         5
7  9.527947         1
8  9.438607         8
9  8.839745         6
10 6.257974         2

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

ここでは話を簡単にするため, R を「直近の購買からの経過時間」でなく,「直近の購買時刻」とします.

モデル化

RF 指標に対して, 非常にシンプルな仮定を置きます.

仮定 1:お客さんの購買行動はレート λポアソンプロセスに従う.
仮定 2:観測は区間 [0, T] で行われた.

この仮定より, 以下が言えます.

仮定 1 よりお客さんの購買行動の間隔は指数分布に従う.
独立に同一の指数分布に従う確率変数の和の分布はガンマ分布になる.

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

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

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

推定

さて, いま興味のある量はポアソンプロセスのレート λ です. λ は顧客の単位時間あたりの購買頻度と解釈できるため, 大きい方がいい数字です.

(1) 式はそのまま尤度と解釈できるので, これを最大化するすることで, λ の推定量 \hat \lambda を得ます.

まず (1) 式の対数をとり,

 \displaystyle \sum _{i=1}^{n} \left\{ x_i \log \lambda + (x_i-1) \log t_i - \log \Gamma (x) - \lambda t_i \right\}

λ微分して

 \displaystyle  \frac{ \sum _{i=1}^{n} x_i}{\lambda} -\sum_{i=1}^{n} t_i

これを0と置いて解くと,

 \displaystyle \hat \lambda =\frac{ \sum_{i=1}^{n} x_i }{\sum_{i=1}^{n} t_i}.

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

Tim <-10 #時刻10まで観測
simu_rf <- function(n=10){
  frequency <- recency <- numeric(n)
  for(i in 1:n){
    N <- rpois(1,0.5*Tim)
    t1 <-sort(runif(N,0,Tim))
    frequency[i] <- N
    recency[i] <- ifelse(N==0,0,t1[N])
  }
  data.frame(recency,frequency)
}
lambda_hat <- numeric(1000)
for(i in 1:1000){
  dat <-simu_rf()
  lambda_hat[i] <- with(dat,sum(frequency)/sum(recency))
}
(mean_lambda_hat <- mean(lambda_hat))
# 0.622057
hist(lambda_hat)
abline(v=1/2,col="red3",lwd=2)
abline(v=mean_lambda_hat,col="blue3",lwd=2)

f:id:abrahamcow:20150831125018p:plain
図は \hat \lambda の値のヒストグラムです. 青い線が推定値の平均, 赤い線が真値です.

推定値にははっきりとバイアスが見て取れます.

なぜでしょうか.

実は (1) 式は完全な尤度になっていませんでした.

直近の購買が時刻 t だった場合, リセンシーは「時刻 t から観測終了時刻 T までの間に, 購買行動が起こっていない」という情報を持っています.

つまりポアソン分布の確率関数より

 \displaystyle \exp(-\lambda (T-t))

を (1) にかけなければいけません.

あらためて尤度は,

 \displaystyle \prod_{i=1}^{n} \frac{ \lambda ^x_i t^{x-1}}{\Gamma(x)} \exp(-\lambda t) \exp(-\lambda (T-t)).

対数をとり, λ で微分すると,

 \displaystyle \sum _{i=1}^{n} \left\{ x_i \log \lambda + (x_i-1) \log t_i - \log \Gamma (x) - \lambda T \right\}

 \displaystyle \sum _{i=1}^{n} \left\{ \frac{x_i}{\lambda}  -  T \right\}

これを0と置いて解くと,

 \displaystyle \tilde \lambda =\frac{ \sum_{i=1}^{n} x_i }{nT}.

この \tilde \lambda についてもシミュレーションを行ってみます.

lambda_tilde <- numeric(1000)
for(i in 1:1000){
  dat <-simu_rf()
  lambda_tilde[i] <- with(dat,sum(frequency)/(10*Tim))
}
(mean_lambda_tilde <- mean(lambda_tilde))
#0.50142
hist(lambda_tilde)
abline(v=1/2,col="red3",lwd=2)
abline(v=mean_lambda_tilde,col="blue3",lwd=2)

f:id:abrahamcow:20150831125815p:plain

パラメータはうまく求まるようです.

 \hat \lambda =\frac{ \sum_{i=1}^{n} x_i }{\sum_{i=1}^{n} t_i} よりも,  \tilde \lambda =\frac{ \sum_{i=1}^{n} x_i }{nT} がいい, というのがちょっと直感に反する感じがして, 個人的にはおもしろかったです.