BTYDplus (version 0.4.0)

pnbd.mcmc.DrawParameters: Hierarchical Bayes variant of Pareto/NBD

Description

pnbd.mcmc.DrawParameters samples parameters via MCMC for a given CBS matrix

Usage

pnbd.mcmc.DrawParameters(cal.cbs, mcmc = 2500, burnin = 500, thin = 50,
  chains = 2, mc.cores = NULL, use_data_augmentation = TRUE,
  param_init = NULL, trace = 100)

Arguments

cal.cbs

data.frame with columns x, t.x, T.cal

mcmc

number of MCMC steps

burnin

number of initial MCMC steps which are discarded

thin

only every thin-th MCMC step will be returned

chains

number of MCMC chains to be run

mc.cores

number of cores to use in parallel (Unix only); defaults to min(chains, detectCores())

use_data_augmentation

determines MCMC method to be used

param_init

list of 2nd-level parameter start values

trace

print logging step every trace iteration

Value

2-element list:

• level_1list of mcmc.list objects; one for each customer, containing individual-level draws

• level_2mcmc.list object containing draws of heterogeneity parameters

Details

method 1) If use_data_augmentation==TRUE MCMC scheme takes advantage of conjugate priors for drawing lambda and mu, by augmentating the parameter space with unobserved lifetime 'tau' and activity status 'z'. See technical appendix to (Abe 2009).

method 2) If use_data_augmentation==FALSE then implementation follows Shao-Hui Ma & Jin-Lan Liu paper http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4344404, i.e. no data augmentation and draws on individual level need to be done via slice sampling. As such it is 10x slower than method 1)

Estimating parameters via Pareto/NBD MCMC can be 10x slower than Pareto/NBD MLE, which itself can be 10x slower than BG/NBD. Both methods exhibit highly autocorrelated draws of r, alpha, s, beta and hence need to be run long, to generate 'enough' draws

References

Ma, Shao-Hui, and Jin-Lan Liu. "The MCMC approach for solving the Pareto/NBD model and possible extensions." Natural Computation, 2007. ICNC 2007. Third International Conference on. Vol. 2. IEEE, 2007. http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4344404

Abe, Makoto. "Counting your customers one by one: A hierarchical Bayes extension to the Pareto/NBD model." Marketing Science 28.3 (2009): 541-553.

See Also

pggg.GenerateData mcmc.DrawFutureTransactions mcmc.PAlive

Examples

# NOT RUN {
set.seed(1)

# generate artificial Pareto/GGG data 
n      <- 1000 # no. of customers
T.cal  <- 32   # length of calibration period
T.star <- 32   # length of hold-out period
params <- list(t=4.5, gamma=1.5, # regularity parameter k_i ~ Gamma(t, gamma)
               r=5, alpha=10,    # purchase frequency lambda_i ~ Gamma(r, alpha)
               s=0.8, beta=12)   # dropout probability mu_i ~ Gamma(s, beta)

data <- pggg.GenerateData(n, T.cal, T.star, params, return.elog=TRUE)
cbs <- data$cbs
elog <- data$elog

# estimate Pareto/NBD MCMC
pnbd.draws <- pnbd.mcmc.DrawParameters(cbs, mcmc=1500, burnin=500, chains=2, thin=50)
plot(pnbd.draws$level_2)

# estimate Pareto/GGG MCMC
pggg.draws <- pggg.mcmc.DrawParameters(cbs, mcmc=1500, burnin=500, chains=2, thin=50)
plot(pggg.draws$level_2, density=FALSE)
rbind("actual"=params, "estimated"=round(summary(pggg.draws$level_2)$quantiles[, "50%"], 2))
#           t   gamma r    alpha s   beta
# actual    4.5  1.5  5    10    0.8  12  
# estimated 3.7  1.2  4.85 9.72  0.71 9.2

coda::gelman.diag(pggg.draws$level_2)
# -> MCMC chains have not converged yet

pggg.mcmc.plotRegularityRateHeterogeneity(pggg.draws)

round(effectiveSize(pggg.draws$level_2))
# -> effective sample size are small for such a short chain

# draw future transaction
pnbd.xstar <- mcmc.DrawFutureTransactions(cbs, pnbd.draws, T.star=cbs$T.star)
pggg.xstar <- mcmc.DrawFutureTransactions(cbs, pggg.draws, T.star=cbs$T.star)

# calculate mean over future transaction draws for each customer
cbs$pnbd.mcmc <- apply(pnbd.xstar, 2, mean)
cbs$pggg.mcmc <- apply(pggg.xstar, 2, mean)

MAPE <- function(a, f) { return(sum(abs(a-f)/sum(a))) }
RMSE <- function(a, f) { return(sqrt(mean((a-f)^2))) }
MSLE <- function(a, f) { return(mean(((log(a+1) - log(f+1)))^2)) }
BIAS <- function(a, f) { return(mean(a)-mean(f)) }
bench <- function(cbs, models) {
  acc <- t(sapply(models, function(model) c(MAPE(cbs$x.star, cbs[, model]),
                                            RMSE(cbs$x.star, cbs[, model]),
                                            MSLE(cbs$x.star, cbs[, model]),
                                            BIAS(cbs$x.star, cbs[, model]))))
  colnames(acc) <- c("MAPE", "RMSE", "MSLE", "BIAS")
  round(acc, 3)
}

bench(cbs, c("pnbd.mcmc", "pggg.mcmc"))
#            MAPE  RMSE  MSLE   BIAS
# pnbd.mcmc  0.532 4.485 0.477 -0.21
# pggg.mcmc 0.488 4.437 0.379 -0.02

# calculate P(active)
cbs$pactive.pnbd.mcmc <- apply(pnbd.xstar, 2, function(x) mean(x>0))
cbs$pactive.pggg.mcmc <- apply(pggg.xstar, 2, function(x) mean(x>0))

# Brier score
c("pnbd.mcmc"=mean((cbs$pactive.pnbd.mcmc-(cbs$x.star>0))^2),
  "pggg.mcmc"=mean((cbs$pactive.pggg.mcmc-(cbs$x.star>0))^2))
#  pnbd.mcmc pggg.mcmc 
# 0.04344944 0.03808333

# calculate P(alive)
cbs$palive.pnbd.mcmc <- mcmc.PAlive(cbs, pnbd.draws)
cbs$palive.pggg.mcmc <- mcmc.PAlive(cbs, pggg.draws)
# }