廿TT

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

阿部誠(2008)RF 指標から生存確率を求める(rstan 版)

モデルと尤度関数

阿部(2004)RF指標から生存確率を求める - 廿TT の続きです。

阿部(2004)によれば, 以下の 2 つの仮定を置くことで, RF(リセンシーとフリクエンシー)指標のみから, 顧客の生存期間を求めることができたのでした.

フリクエンシーを x, リセンシー(直近購買日)を t で表すことにします。

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

潜在変数 z は顧客が時刻 T の時点で生存しているときに 1, 離脱しているときに 0 の値を取るとします.

y は顧客が離脱した時点とします. これも観測されません.

z=1 のとき, 対数尤度は,

x\log(\lambda)-(\lambda+\mu)T

z=0 のとき, 対数尤度は,

x\log(\lambda)+\log(\mu)-(\lambda+\mu)y

となります.

z が値 1 を取る確率は,

\displaystyle p = \frac{1}{1+\frac{\mu}{\lambda+\mu}( e^{(\lambda+\mu)(T-t)}) }

です.

詳しくは論文をみてください.

Stan

阿部(2008)では λμ の対数に(独立でない)二変量正規分布を仮定していたりしますが、ここでは一様分布にします.

y区間 [t, T] の切断指数分布を仮定します.

問題は z ですがこれは離散パラメータなので周辺化します.

取りうるすべての値 (0, 1) を足し合わせて消去します.

RStanで離散パラメータを含むモデルの推定(disaster model) - 廿TT などを参照してください.

Stan コードは以下のようになりました.

data {
  real<lower=0> T;
  int<lower=0> x;
  real<lower=0> t; 
}
parameters {
  real<lower=0> lambda;
  real<lower=0> mu;
  real<lower=t,upper=T> y;
}
transformed parameters{
  real<lower=0,upper=1> p;
  p = 1/(1+mu*( exp((lambda+mu)*(T-t))/(lambda+mu) ));
}
model {
  target += bernoulli_lpmf(1|p)+ #z=1 のとき
    x*log(lambda)-(lambda+mu)*T;
  target += bernoulli_lpmf(0|p)+ #z=0 のとき
    x*log(lambda)+log(mu)-(lambda+mu)*y;
  y~exponential(mu)T[t,T];
}

シミュレーションで生成したデータ・セットから、パラメータを推定してみます.

f:id:abrahamcow:20160806001233p:plain

トレースラインをみると定常分布に収束しているようです.

次のヒストグラムはパラメータ λ の事後分布です. 最尤法による推定値を青い線, 真値をオレンジ, 事後平均を赤い線で表しています.

f:id:abrahamcow:20160806001306p:plain

次のヒストグラムはパラメータ μ の事後分布です.

f:id:abrahamcow:20160806001309p:plain

これを見る限りでは最尤法のほうが真値に近づくようです.

以下に R のコードを貼っておきます.

シミュレーションと最尤推定については, 阿部(2004)RF指標から生存確率を求める - 廿TT を参照してください.

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,]
}
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
abemodel1 <-stan_model("~/Documents/abe2008.stan")

dat0 <-simu_rf2(10) #多めに出してる
datlist<- list(t=dat0$recency[1]-dat0$first[1],
               T=Tim-dat0$first[1], x= dat0$frequency[1])

fit <- sampling(abemodel1,datlist)
traceplot(fit)
get_posterior_mean(fit)

dat <- dat0[1,]
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 <-optim(c(1,1),objfun,control=list(fnscale=-1))
lambdaMLE <- exp(res$par[1])
muMLE <- exp(res$par[2])
ex <- extract(fit)
hist(ex$lambda)
abline(v=c(lambdaMLE,mean(ex$lambda),lambda),col=c("blue","red","orange2"),lwd=2)
hist(ex$mu)
abline(v=c(muMLE,mean(ex$mu),mu),col=c("blue","red","orange2"),lwd=2)

追記

と書きましたがすいません Stan コードがまちがっていたようです.

尤度を足し合わせて周辺化するので、対数尤度の和ではなく log_sum_exp を取らないといけなかった.

シミュレーションをやり直しました.

data {
  real<lower=0> T;
  int<lower=0> x;
  real<lower=0> t; 
}
parameters {
  real<lower=0> lambda;
  real<lower=0> mu;
  real<lower=t,upper=T> y;
}
transformed parameters{
  real<lower=0,upper=1> p;
  p = 1/(1+mu*( exp((lambda+mu)*(T-t))/(lambda+mu) ));
}
model {
  real lp[2];
  lp[1] = bernoulli_lpmf(1|p)+
    x*log(lambda)-(lambda+mu)*T;
  lp[2] = bernoulli_lpmf(0|p)+
    x*log(lambda)+log(mu)-(lambda+mu)*y;
  target+=log_sum_exp(lp);
  y~exponential(mu)T[t,T];
}

f:id:abrahamcow:20160806085659p:plain

f:id:abrahamcow:20160806085713p:plain

f:id:abrahamcow:20160806085717p:plain