RとStanでPareto/NBD, Gamma-Gammaモデル

はじめに

購買データの分析と言えば、RFM分析が有名ですね。リーセンシ―Recency（直近購買からの経過時間）、フリークエンシーFrequency（購買回数）、マネタリーMonetary-value（平均購買額)は購買行動を簡潔に表す指標です。これらの指標と確率分布を利用して個人レベルの購買行動や特性を予測・分析するモデルがマーケティングでは様々考案されています。このアプローチは、顧客の生涯価値（Customer Lifetime Value, CLV）を考え、個人レベルの購買力総合指標を用いた顧客管理に非常に役立ちます。

RFMではダメなのか？

RFM自体は顧客の潜在的な（真の）購買特性を表すものではなく、特性から発生した行動の間接的な結果です。それゆえ顧客の把握という意味では欠点があります。例えば同じリーセンシーのAさん、Bさんがいたとしてフリークエンシーが異なれば、AさんとBさんの価値は同じにはなりません。RFMには複雑な相互関係が発生しています。
そこで、潜在的な購買特性としてPLS（購買頻度：Purchase rate、生存期間：Lifetime、1回当たりの購買金額：Spending per transaction）があります。相互関係を完全には補正できなくとも、RFMの欠点を補っているといえます。よってPLSをモデルから求め、CLVを算出してやります。

モデル概要

RFMデータ

顧客はサービスに登録してから離脱するまでの期間に購買が行われています。

顧客の購買履歴データの表記　阿部(2007)より
RFMデータで表すと、Recency = T-t、frequency = x + 1, Monetaryは各xでの額を平均した値ですね。

モデル

さて、上記したPLSを導くために、この期間中Tに行われた購買行動をモデル化します。モデルは”購買の発生”と”購買額”の2つの別々のモデルで成り立っています。今回行うのは、Fader, Hardie and Lee(2005)にて提唱された ”購買の発生”をPareto/NBD, ”購買額”をGamma-Gammaと呼ばれる仮定がもとになっているモデルです。
またそのモデル（とPLS)は、各顧客のRFMデータを用いた3つの単純な購買パラメータで表現できます。

• λ：購買頻度(Transaction_rate)・・・顧客毎の平均購買回数（P）
• μ：離脱率(Dropout_rate)・・・顧客毎にサービスから離脱していなくなる確率、生存期間（L）を出すのに必要
• as：平均購買額(Average_spending)・・・顧客毎の平均購買金額（S）

個人の異質性（差）

このモデルは個人の行動だけでなく、個人の差（各パラメータの違い）も確率分布で表すところにあります。
つまり、

• RFMデータがあるパラメータ（λ, μ等）を持つ確率分布から発生している
• 各パラメータの個人の違いも確率分布で表している

購買モデル（Pareto/NBD）

購買（transaction）の発生を表すモデルです。

仮定①：顧客が生存中という条件の元、特定の期間（初期購買~t）中のリピート購買回数xが購買頻度λのポアソン過程に従って決まる
仮定②：顧客の生存期間（初期購買～離脱の期間：τ）は離脱率μの指数分布に従って決まる

• ①：P(X(t) = x) ~ Poisson(λt)
• ②：τ ~ Exponential(μe^(-μτ))

購買額モデル（gamma-gamma）

平均購買額(Average_spending)を表すモデルです。

仮定③：各顧客の毎回の購買額z={z1, z2, ... zi}はガンマ分布に従って決まる

• ③：zi ~ Gamma(p, ν)　➡　as ~ Gamma(px, νx)
• v ~ Gamma(q, y)

as = \sum_{i=1}^{x} z_i/x 　の平均購買額パラメータ

個人の異質性（差）モデル

各モデルのパラメータの個人差を表す仮定です。
個人差仮定：個人ごとに異なるパラメータλとμがガンマ分布に従う

• λ ~ Gamma(rr, α)
• μ ~ Gamma(ss, β)

CLV

λ, μ, zを利用して、購買回数と購買額の期待値や総合指標であるCLVを計算します。特定期間のCLVを単純に求めても良いかもしれません。

E[x] = \frac{λ}{μ} * (1 - e^{−μt}) ・・・期待リピート購買回数\\\\
CLV = λ * E[z] / (μ + δ) ・・・顧客生涯価値。δ = (1 + d)\\\\
 \\\\

termLTV = E[x] * E[z]・・・特定の指定期間中の期待購買額

dは年間割引率で購買力の経過時間減少を表しています。

RStanで実行

パラメータの推定と言えば、Stanですね。行動-個人差のような階層性のあるモデルはベイズ向きと言えます（もちろんRのライブラリでも確率分布のパラメータ推定は可能です）。しかし一番の理由は個人のCLVを算出するのに必要なλ、μを個人レベルで推定できることです。
データは本モデルの論文でも使われているCDNOWサンプルデータセットです。CDNOWのウエブサイトで買われたCDの顧客別購買履歴78週分（18ヶ月）です。Fader, Hardie and Leeと阿部同様最初の39週でモデルを作成（推定）し、最後の39週分で検証を行います。Stanはサンプル数が多いと実行時間がかかるので、減らして実行しています。

stan　nbd_gamma_gamma01

//Qiitaのコードブロック、Stanは未対応ですね。
data {
  int<lower = 1> N;            // サンプル数：N
  int<lower = 0> Term_m;       // テストに使用する期間
  vector<lower = 0>[N] X;      // 再購買回数：Repeat transactions (frequency)
  vector<lower = 0>[N] Tx;     // 0 ~ t
  vector<lower = 0>[N] Ti;     // term : 0 ~ T
  vector[N] AS;                // 平均購買額：Average Spending
}

transformed data {
  vector<lower = 0>[N] X_total;  // 総購買回数
  X_total = X + 1;
}

parameters {
  vector<lower = 0>[N] lambda; // 購買頻度率：Transaction rate
  vector<lower = 0>[N] mu;     // 離脱確率：Dropout probability
  // lambda gamma分布のパラメータ：購買
  real<lower = 0> rr;          // 購買λのshape：Transaction shape parameter
  real<lower = 0> alpha;       // 購買λのscale：Transaction scale parameter
  // mu gamma分布のパラメータ：離脱
  real<lower = 0> ss;          // 離脱率μのshape：Dropout shape parameter
  real<lower = 0> beta;        // 離脱率μのscale：Dropout scale parameter
  // as gamma分布のパラメータ：額
  real <lower=0> p;            // 額のshape gamma        
  vector<lower=0>[N] v;        // 額のscale gamma
  // 事前分布のパラメータ
  real <lower=0> q;            // 額のscale vのshape gamma  
  real <lower=0> y;            // 額のscale vのscale gamma    
}

transformed parameters {
  vector[N] log_lambda; 
  vector[N] log_mu;
  vector[N] log_mu_lambda;
  vector<lower=0> [N] px;   
  vector<lower=0> [N] nx;  
  vector[N] p_1;
  vector[N] p_2;
  // repeat購買と離脱率 repeat parameter
  log_lambda = log(lambda);
  log_mu = log(mu);
  log_mu_lambda = log(mu + lambda);
  // 支出額のgamma distribution parameter
  px = p * X_total;  
  nx = v .* X_total;
  // L(λ, μ | X, Tx, T) 
  p_1 = X .* log_lambda + log_mu - log_mu_lambda - Tx .* mu - Tx .* lambda; 
  p_2 = X .* log_lambda + log_lambda - log_mu_lambda - Ti .* mu - Ti .* lambda; 

}

model {
  // repeat購買と離脱率の分布と推定
  lambda ~ gamma(rr, alpha);
  mu ~ gamma(ss, beta);
  rr ~ exponential(1);
  alpha ~ exponential(1);
  ss ~ exponential(0.1);
  beta ~ exponential(0.1);
  target += log(exp(p_1) + exp(p_2));
  // 額の事前分布 Priors for spend
  p ~ exponential(0.1);    
  q ~ exponential(0.1);    
  y ~ exponential(0.1);    
  v ~ gamma(q, y); 
  // 額の推定 Likelihood for spend
  AS ~ gamma(px, nx); 
}

generated quantities {
  vector[N] p_alive;        // T時点の生存確率：Probability that they are still "alive"    
  vector[N] exp_trans;      // Expected number of transactions 
  vector[N] as_pred;        // Per transaction spend
  vector[N] ltv;         // LTV Lifetime value
  vector[N] clv;         // CLV customer Lifetime value
  vector[N] exp_transM1;
  vector[N] exp_transM2;
  vector[N] exp_transM3;
  vector[N] exp_transM4;
  vector[N] exp_transM5;
  vector[N] exp_transM6;
  vector[N] exp_transM7;
  vector[N] exp_transM8;
  vector[N] exp_transM9;
  vector[N] exp_spendM1;
  vector[N] exp_spendM2;
  vector[N] exp_spendM3;
  vector[N] exp_spendM4;
  vector[N] exp_spendM5;
  vector[N] exp_spendM6;
  vector[N] exp_spendM7;
  vector[N] exp_spendM8;
  vector[N] exp_spendM9;
  
  for(i in 1:N) {
    p_alive[i] = 1/(1+mu[i]/(mu[i]+lambda[i])*(exp((lambda[i]+mu[i])*(Ti[i]-Tx[i]))-1)); //T時点の将来の期待生存確率 P(τ>T)
    exp_trans[i] = (lambda[i]/mu[i])*(1 - exp(-mu[i]*Term_m)); //将来の期待購買回数 E[x]
    as_pred[i] = gamma_rng(px[i], nx[i]);  //将来の1回当たりの平均期待購買金額 E[Sn]
    
    exp_transM1[i] = (lambda[i] / mu[i]) * (1 - exp(-mu[i] * 1)); 
    exp_transM2[i] = (lambda[i] / mu[i]) * (1 - exp(-mu[i] * 2)); 
    exp_transM3[i] = (lambda[i] / mu[i]) * (1 - exp(-mu[i] * 3)); 
    exp_transM4[i] = (lambda[i] / mu[i]) * (1 - exp(-mu[i] * 4)); 
    exp_transM5[i] = (lambda[i] / mu[i]) * (1 - exp(-mu[i] * 5)); 
    exp_transM6[i] = (lambda[i] / mu[i]) * (1 - exp(-mu[i] * 6)); 
    exp_transM7[i] = (lambda[i] / mu[i]) * (1 - exp(-mu[i] * 7)); 
    exp_transM8[i] = (lambda[i] / mu[i]) * (1 - exp(-mu[i] * 8)); 
    exp_transM9[i] = (lambda[i] / mu[i]) * (1 - exp(-mu[i] * 9));
    
    exp_spendM1[i] = exp_transM1[i] * as_pred[i];
    exp_spendM2[i] = exp_transM2[i] * as_pred[i];
    exp_spendM3[i] = exp_transM3[i] * as_pred[i];
    exp_spendM4[i] = exp_transM4[i] * as_pred[i];
    exp_spendM5[i] = exp_transM5[i] * as_pred[i];
    exp_spendM6[i] = exp_transM6[i] * as_pred[i];
    exp_spendM7[i] = exp_transM7[i] * as_pred[i];
    exp_spendM8[i] = exp_transM8[i] * as_pred[i];
    exp_spendM9[i] = exp_transM9[i] * as_pred[i];
  }
  clv = lambda .* as_pred ./ (mu + 0.0124) ; //年間割引率15%、月単位
  term_ltv = exp_trans .* as_pred;
}

}

使用するデータの用意と実行はRです。

R

library(tidyverse)
cdnow_sample <- read.table(file=paste0('CDNOW_sample/CDNOW_sample.txt'))

week39 = as.Date('1997-01-01') +7*39
custs <- cdnow_sample %>%
  select(-V1) %>%
  mutate(customer_id = V2,
         logged = as.Date(as.character(V3), format = "%Y%m%d"),
         total = V5) %>%
  group_by(customer_id, logged) %>%
  summarise(spending = sum(total)) %>% ungroup() %>%
  arrange(customer_id,logged)
custs$week <- ceiling( (as.numeric(date(custs$logged)) -9861) /7)
custs$month <- ceiling( (as.numeric(date(custs$logged)) -9861) /30.41)

start_day <- date(min(custs$logged))
end_day <- date('1997-10-01') # 39week

customer_data <- custs %>%
  filter(logged < '1997-10-01') %>%
  group_by(customer_id) %>%
  summarise(x = n() - 1,
            sday = start_day,
            mind =min(logged),
            maxd =max(logged),
            first_pay = first(spending),
            t_x = as.numeric(difftime(max(logged), min(logged), units = "days")/30.41),
            t_i = as.numeric(difftime(end_day, min(logged), units = "days")/30.41),
            as = mean(spending)) %>% ungroup() %>%
  mutate(Ti=t_i-t_x, ncustomer_id = row_number() ) %>%
  arrange(ncustomer_id)

customer_data100 <- customer_data %>% filter(ncustomer_id <= 100)

model_data <- list(
  X = customer_data100$x,
  Tx = as.numeric(customer_data100$t_x),
  Ti = as.numeric(customer_data100$t_i),
  AS = customer_data100$as,
  N = nrow(customer_data100),
  Term_m = 9
)


library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

stanmodel <- stan_model(file='nbd_gamma_gamma01.stan')

fit<-sampling(stanmodel,data=model_data, seed=1234,
              chains=4, iter=3000,warmup=1000,thin=1
)
# コードの確認
get_stanmodel(fit)
# サンプル抽出
ms <- rstan::extract(fit)
# 実行サンプル数 = (iter - warmup) × chain / thin)
N_mcmc <- length(ms$lp__)

summary_fit<-data.frame(summary(fit)$summary[,c("mean","se_mean","sd","2.5%","25%","50%","75%","97.5%","n_eff","Rhat")]) #Rhat<1.1以下で収束
summary_fit<-summary_fit %>% mutate(parameters= rownames(summary_fit))
Rhat1.1<-summary_fit%>%filter(Rhat>1.1| is.na(Rhat))

結果と検証

①期間-累積購買回数

検証期間中の顧客の累積購買回数を時系列に実績と予測で比較します。

予測の方が少々購買回数を多く見積もる。

②検証期間中の購買回数ヒストグラム

検証期間中の購買回数別顧客数（ヒストグラム）を実績と予測で比較します。

0回実績が大半なので予測も大半が0回になれば精度が必然的に上がってしまう。課題としては0回ではない顧客数と購買回数を予測することが重要と思える。

③購買回数セグメント別平均購買回数

モデル推定期間中の購買回数で顧客をセグメント分けし、検証期間中の購買回数の実績と予測で比較します。

お、いいかんじ

④購買回数セグメント別平均購買額

モデル推定期間中の購買回数で顧客をセグメント分けし、検証期間中の購買額の実績と予測で比較します。

こちらも、いいかんじ

とはいえ

検証期間中の購買が0回の顧客が多いため、どうしてもそれに引っ張られてしまいますね。顧客管理を目的として顧客の購買力を把握するのであれば、モデル推定期間中の購買が複数回の顧客に絞ってモデルを検討してみる必要がありそうです。

Rでパラメータ推定する（非線形方程式を解く方法）

Stanはサンプル数が増えるととにかく時間がかかるようになります。個人レベルのパラメータであるλとμがサクッと推定できればCLVもすぐ計算できます。
BTYDライブラリでは個別のパラメータ推定はできないので、別の方法で推定してStanの結果と比較してみました。

R

# --------------------------------------------------------------------
# λとμの推定（Rで）
# --------------------------------------------------------------------
loglik <- function(param, Ti,Tx,X){
  p_1 = X * log(param[1]) + log(param[2]) - log(param[2]+param[1]) - Tx * param[2] - Tx * param[1]; #生存：P(λ, μ|τ>T)
  p_2 = X * log(param[1]) + log(param[1]) - log(param[2]+param[1]) - Ti * param[2] - Ti * param[1]; #死亡：P(λ, μ|Tx<τ<T)
  return(p_1+p_2)
}

N<-length(model_data$X)
params <- c()
for(i in 1:N){

  Ti <- customer_data100$t_i[i]
  Tx <- customer_data100$t_x[i]
  X <- customer_data100$x[i]
  opt <- optim(c(0.1,0.1), loglik, Ti=Ti, Tx=Tx, X=X,control=list(fnscale=-1),method="BFGS")
  params <- rbind(params,c(opt$par,opt$convergence))
  # opt <- c()
}

# --------------------------------------------------------------------
# 結果比較
# --------------------------------------------------------------------
# ----------------------------------
# Stanのlambda mu
# ----------------------------------
lambda <- summary_fit %>% 
  filter(str_detect(parameters,'lambda'))  %>%
  filter(str_detect(parameters,'log_lambda', negate = TRUE))  %>%
  filter(str_detect(parameters,'log_mu_lambda', negate = TRUE))  %>%
  mutate(customer_id=customer_data100$customer_id) %>%
  group_by(customer_id) %>%
  summarise(stan_med_lambda =first(X50.))

mu <- summary_fit %>% 
  filter(str_detect(parameters,'mu'))  %>%
  filter(str_detect(parameters,'log_mu', negate = TRUE))  %>%
  filter(str_detect(parameters,'log_mu_lambda', negate = TRUE))  %>%
  mutate(customer_id=customer_data100$customer_id) %>%
  group_by(customer_id) %>%
  summarise(stan_med_mu =first(X50.))

lambda_mu <-lambda %>% full_join(mu,by='customer_id')

lambda_mu$r_lambda <- params[,1]
lambda_mu$r_mu <- params[,2]

ggplot(data= lambda_mu) +
  aes(x=stan_med_lambda,y=r_lambda) +
  geom_point() +
  xlim(0,1.5) +
  ylim(0,1.5)

ggplot(data= lambda_mu) +
  aes(x=stan_med_mu,y=r_mu) +
  geom_point() +
  xlim(0,0.1) +
  ylim(0,0.1)

StanとRの比較

stanの50%値とRの推定結果が同じなら45度の直線関係になるはず。

λの比較

割ときれいに直線。

μの比較

μはかなりずれてますね・・・下の方でたまっているサンプルはリピート購買回数Xが0回の顧客です。やはりXの回数不足が否めません。

おわりに

顧客の購買を説明する面白い考え方だと思います。モデルの詳しい説明や考え方は阿部(2007)が参考になります。とはいえ、現場レベルでどう使うのか？という部分ではPLSの考え方を利用したやずやさんのマーケティング手法が楽だったりしますね。