Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
/* -----------------------------------------------------------------------------
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)) ;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment