rfm_hierarchical1702.stan

/* -----------------------------------------------------------------------------
Estimate latent variable hierarchical Bayes RFM model 

stan ver. 2.14

CREATED: by ill-identified, 18/APR/2016 
REVISED: 25/FEB/2017, REVISED
------------------------------------------------------------------------------- */

data {
  int<lower=1>              N ; // num of customers
  int<lower=1>              K ; // max number of spending for each customers
  vector<lower=0>[N]     Time ; // duration between the first observed purchase and the end of monitoring
  vector<lower=0>[N]     time ; // duration between fisrt and last purchase 
  vector<lower=0>[N]        x ; // frequency (number of repeat)
  int<lower=1>      NSpend[N] ; // number of spending for each customers
  matrix<lower=0>[N,K]  Spend ; // spendings
  real<lower=0, upper=1> delta ; // discount rate for CLV
}
transformed data{
  vector<lower=0>[N] R ; // recency
  R = Time - time ;
}

parameters {
  vector[3]  param3[N] ; // [lambda, mu, eta] as vector in **log**
  // param3[, 1]:lambda, param3[, 2]: mu, param3[, 3]: eta
  vector<lower=0, upper=1>[N] tau_rand ; // intermediate latent indicator
  vector<lower=0>[N]             omega ; // variance of monetary

  // hyperparameter
  vector[3]      theta0 ;         // mean of 3 main params
  cholesky_factor_corr[3] Gamma0 ; // corr of 3 main params
  real<lower=1>  lkj ;        // corr hyperparameter
}

transformed parameters{
  vector[N] tau ;
  vector<lower=0, upper=1>[N] xi ;
  tau = time + R .* tau_rand ;
  // xi = P(zeta = 1), or if surviving after Time[i] )
  for ( i in 1:N){
    if(R[i] == 0){
      xi[i] = 1 ;
    }
    else {
      xi[i] = 1 / fma(exp(param3[i][2])/(exp(param3[i][1]) + exp(param3[i][2])), 
                    expm1((exp(param3[i][1])+exp(param3[i][2])) * R[i]), 1) ;
    }
  }
}

model {
  /* hyper priors */
  lkj ~ student_t(4, 1, 10) T[1, ] ;
  tau_rand ~ uniform(0, 1) ;
  Gamma0 ~ lkj_corr_cholesky(lkj) ;
  theta0 ~ normal(0, 20) ;
  // prior dist.
  param3 ~ multi_normal_cholesky(theta0, Gamma0) ;

  for(i in 1:N)
    omega[i] ~ student_t(4, 0, 10) T[0, ] ;
  // latent duration
  // improper prior
  target += exponential_lpdf(tau | exp(to_vector(param3[, 2])) ) ;

  for (i in 1:N) {
    if(x[i] > 0)
      target += x[i] * param3[i][1] + lmultiply(x[i] - 1, time[i]) - lgamma(x[i]) ;

    // marginalize out zeta
    target += log_sum_exp(
      bernoulli_lpmf(1| xi[i]) - (exp(param3[i][1]) + exp(param3[i][2]) ) * Time[i],
      bernoulli_lpmf(0| xi[i]) + param3[i][2] - (exp(param3[i][1]) + exp(param3[i][2])) * tau[i]
    ) ;

    // spending ~ lognormal
    for( k in 1:NSpend[i])
      Spend[i,k] ~ lognormal(param3[i][3], omega[i]) ;
  }
}


generated quantities{
  vector[N]     lambda ; // frequency parameter
  vector[N]         mu ; // survival parameter
  vector[N]        eta ; // monetary value parameter
  vector[N]        CLV ; // CLV of each customer
  vector[N]        Tau ; // to evaluate distribution

  lambda  = exp(to_vector(param3[,1])) ;
  mu      = exp(to_vector(param3[,2])) ;
  eta     = exp(to_vector(param3[,3])) ;

  CLV = (lambda .* eta .* square(omega)) ./ (2 * (mu + delta)) ;
}