BTYDplus (version 0.4.0)

pggg.mcmc.DrawParameters: Hierarchical Bayes implementation of Pareto/GGG

Description

Returns 2-element list level_1: 3-dim array [draw x parameter x cust] wrapped as coda::mcmc.list object level_2: 2-dim array [draw x parameter] wrapped as coda::mcmc.list object

Usage

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

Arguments

cal.cbs

data.frame with columns x, t.x, T.cal, litt; e.g. output of elog2cbs

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())

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

See Also

link{pggg.GenerateData} mcmc.PAlive mcmc.DrawFutureTransactions elog2cbs

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)
# }