RFM_Bayes.stan

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

stan ver. 2.9

CREATED BY: 09/APR/2015, ill-identified 
------------------------------------------------------------------------------- */

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 
  int<lower=0>           x[N] ; // 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(x)**
  // 1: lambda, 2: mu, 3: eta

  vector<lower=0,upper=1>[N]   xi ; // latent indicator variable 
  vector<lower=0>[N]          tau ; // latent duration variables
  vector<lower=0>[N]        omega ; // variance of monetary

  // hyperparameter
  vector[3]      theta0 ;
  cov_matrix[3]  Gamma0 ;
}

transformed parameters{
  vector[N] zeta ;
  // xi ~ U(0,1), zeta = 1(xi < p(zeta=1) )
  for ( i in 1:N){
    zeta[i] <- 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 parameters */
  vector[N] a;
  Gamma0 ~ inv_wishart(3, 10*diag_matrix(rep_vector(1,3)) ) ;
  theta0 ~ uniform(-10, 10) ;
  // prior dist.
  param3 ~ multi_normal(theta0, Gamma0) ;
  omega ~ uniform(0,10) ;
  xi ~ uniform(0,1) ;
  
  /* --- likelihood --- */
  for ( i in 1:N) {
    if(x[i] > 0)
      increment_log_prob(x[i]*param3[i][1] + multiply_log(x[i]-1, time[i]) - lgamma(x[i]) ) ;
    // spending ~ lognormal
    for( k in 1:NSpend[i])
      Spend[i,k] ~ lognormal(param3[i][3], omega[i]) ;
      
    // zeta-dependending part
    if (zeta[i] == 1 || R[i] == 0 ){ // surviving during monitoring 
      increment_log_prob(-(exp(param3[i][1]) + exp(param3[i][2]) )*Time[i] ) ;
    }
    else if (zeta[i] == 0){
      tau[i] ~ exponential(exp(param3[i][1])+exp(param3[i][2]) ) T[,Time[i]];
      //increment_log_prob( exponential_log(tau[i], exp(param3[i][1])+exp(param3[i][2]) ) ) ;
      /*increment_log_prob(exponential_log(tau[i], exp(param3[i][1])+exp(param3[i][2]) ) 
                         -log(exponential_cdf(Time[i], exp(param3[i][1])+exp(param3[i][2])) - 
                              exponential_cdf(time[i], exp(param3[i][1])+exp(param3[i][2]))) ) ;*/
      increment_log_prob(param3[i][2] - (exp(param3[i][1]) + exp(param3[i][2])) * tau[i]) ;
    }
    else{
      reject("zeta is neither 0 nor 1 !") ;
    }
  }
}

generated quantities{
  vector[N]    log_lik ; // log likelihood ;
  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]      pzeta ;

  lambda  <- exp(to_vector(param3[,1])) ;
  mu      <- exp(to_vector(param3[,2])) ;
  eta     <- exp(to_vector(param3[,3])) ;
  
  for( i in 1:N){
    // P(zeta=1)
    pzeta[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) ;
    // log likelihood 
    log_lik[i] <- x[i]*param3[i][1] + (x[i]-1)*log(time[i])  ;
    if(x[i] > 0)
      log_lik[i] <- log_lik[i] - lgamma(x[i]) ;
    if (zeta[i] == 1)
      log_lik[i] <- log_lik[i] -(exp(param3[i][1])+param3[i][2])*Time[i] ;
    else
      log_lik[i] <- log_lik[i] + param3[i][2] - (exp(param3[i][1]) + exp(param3[i][2])) * tau[i] ;
      
    for(k in 1:NSpend[i])
      log_lik[i] <- log_lik[i] + lognormal_log(Spend[i,k], eta, omega[i]) ;
      
    // CLV 
    CLV[i] <- (lambda[i]*eta[i]*omega[i]^2) / (2*(mu[i] + delta) ) ;
  }
}