Nothing
################################################## BG/NBD estimation, visualization functions
library(hypergeo)
# Two things discovered in this script so far:
# -- bgnbd.cbs.LL should be called with the un-compressed version of cal.cbs, the 3-column one
# -- bgnbd.LL spec, as written, won't avoid the large x problem. Patched that, not tested yet.
#' Define general parameters
#'
#' This is to ensure consistency across all functions that require common bits
#' and bobs.
#'
#' @inheritParams bgnbd.LL
#' @inheritParams bgnbd.ConditionalExpectedTransactions
#' @param func function calling dc.InputCheck
#' @param hardie if TRUE, use \code{\link{h2f1}} instead of
#' \code{\link[hypergeo]{hypergeo}} when you call this function from within
#' \code{\link{bgnbd.ConditionalExpectedTransactions}}.
#' @return a list with things you need for \code{\link{bgnbd.LL}},
#' \code{\link{bgnbd.PAlive}} and
#' \code{\link{bgnbd.ConditionalExpectedTransactions}}
#' @seealso \code{\link{bgnbd.LL}}
#' @seealso \code{\link{bgnbd.PAlive}}
#' @seealso \code{\link{bgnbd.ConditionalExpectedTransactions}}
bgnbd.generalParams <- function(params,
func,
x,
t.x,
T.cal,
T.star = NULL,
hardie = NULL) {
inputs <- try(dc.InputCheck(params = params,
func = func,
printnames = c("r", "alpha", "a", "b"),
x = x,
t.x = t.x,
T.cal = T.cal))
if('try-error' == class(inputs)) return(inputs)
x <- inputs$x
t.x <- inputs$t.x
T.cal <- inputs$T.cal
r <- params[1]
alpha <- params[2]
a <- params[3]
b <- params[4]
# last two components for the alt specification
# to handle large values of x (Solution #2 in
# http://brucehardie.com/notes/027/bgnbd_num_error.pdf,
# LL specification (4) on page 4):
C3 = ((alpha + t.x)/(alpha + T.cal))^(r + x)
C4 = a / (b + x - 1)
# stuff you'll need in sundry places
out <- list()
out$PAlive <- 1/(1 + as.numeric(x > 0) * C4 / C3)
# do these computations only if needed: that is,
# if you call this function from bgnbd.LL
if(func == 'bgnbd.LL') {
# a helper for specifying the log form of the ratio of betas
# in http://brucehardie.com/notes/027/bgnbd_num_error.pdf
lb.ratio = function(a, b, x, y) {
(lgamma(a) + lgamma(b) - lgamma(a + b)) -
(lgamma(x) + lgamma(y) - lgamma(x + y))
}
# First two components -- D1 and D2 -- for the alt spec
# that can handle large values of x (Solution #2 in
# http://brucehardie.com/notes/027/bgnbd_num_error.pdf)
# Here is the D1 term of LL function (4) on page 4:
D1 = lgamma(r + x) -
lgamma(r) +
lgamma(a + b) +
lgamma(b + x) -
lgamma(b) -
lgamma(a + b + x)
D2 = r * log(alpha) - (r + x) * log(alpha + t.x)
# original implementation of the log likelihood
# A = D2 + lgamma(r + x) - lgamma(r)
# B = exp(lb.ratio(a, b + x, a, b)) *
# C3 +
# as.numeric((x > 0)) *
# exp(lb.ratio(a + 1, b + x - 1, a, b))
# out$LL = sum(A + log(B))
# with the corection for avoiding the NUM! problem:
out$LL = D1 + D2 + log(C3 + as.numeric((x > 0)) * C4)
}
# if T.star is not null, then this can produce
# conditional expected transactions too. this is
# another way of saying that you are calling this
# function from bgnbd.ConditionalExpectedTransactions,
# in which case you also need to set hardie to TRUE or FALSE
if(!is.null(T.star)) {
stopifnot(hardie %in% c(TRUE, FALSE))
term1 <- (a + b + x - 1) / (a - 1)
if(hardie == TRUE) {
hyper <- h2f1(r + x,
b + x,
a + b + x - 1,
T.star/(alpha + T.cal + T.star))
} else {
hyper <- Re(hypergeo(r + x,
b + x,
a + b + x - 1,
T.star/(alpha + T.cal + T.star)))
}
term2 <- 1 -
((alpha + T.cal)/(alpha + T.cal + T.star))^(r + x) *
hyper
out$CET <- term1 * term2 * out$PAlive
}
out
}
#' BG/NBD Log-Likelihood
#'
#' Calculates the log-likelihood of the BG/NBD model.
#'
#' \code{x}, \code{t.x} and \code{T.cal} may be vectors. The standard rules for
#' vector operations apply - if they are not of the same length, shorter vectors
#' will be recycled (start over at the first element) until they are as long as
#' the longest vector. It is advisable to keep vectors to the same length and to
#' use single values for parameters that are to be the same for all
#' calculations. If one of these parameters has a length greater than one, the
#' output will be also be a vector.
#'
#' @param params BG/NBD parameters - a vector with r, alpha, a, and b, in that
#' order. r and alpha are unobserved parameters for the NBD transaction
#' process. a and b are unobserved parameters for the Beta geometric dropout
#' process.
#' @param x number of repeat transactions in the calibration period T.cal, or a
#' vector of transaction frequencies.
#' @param t.x time of most recent repeat transaction, or a vector of recencies.
#' @param T.cal length of calibration period, or a vector of calibration period
#' lengths.
#'
#' @seealso \code{\link{bgnbd.EstimateParameters}}
#' @seealso \code{\link{bgnbd.cbs.LL}}
#'
#' @return A vector of log-likelihoods as long as the longest input vector (x,
#' t.x, or T.cal).
#'
#' @examples
#' data(cdnowSummary)
#'
#' cal.cbs <- cdnowSummary$cbs
#' # cal.cbs already has column names required by method
#'
#' # random assignment of parameters
#' params <- c(0.5, 6, 1.2, 3.3)
#' # returns the log-likelihood of the given parameters
#' bgnbd.cbs.LL (params, cal.cbs)
#'
#' # compare the speed and results to the following:
#' cal.cbs.compressed <- dc.compress.cbs(cal.cbs)
#' bgnbd.cbs.LL(params, cal.cbs.compressed)
#'
#' # Returns the log likelihood of the parameters for a customer who
#' # made 3 transactions in a calibration period that ended at t=6,
#' # with the last transaction occurring at t=4.
#' bgnbd.LL(params, x=3, t.x=4, T.cal=6)
#'
#' # We can also give vectors as function parameters:
#' set.seed(7)
#' x <- sample(1:4, 10, replace = TRUE)
#' t.x <- sample(1:4, 10, replace = TRUE)
#' T.cal <- rep(4, 10)
#' bgnbd.LL(params, x, t.x, T.cal)
bgnbd.LL <- function(params,
x,
t.x,
T.cal) {
bgnbd.generalParams(params = params,
func = 'bgnbd.LL',
x = x,
t.x = t.x,
T.cal = T.cal)$LL
}
#' BG/NBD Log-Likelihood Wrapper
#'
#' Calculates the log-likelihood sum of the BG/NBD model.
#'
#' Note: do not use a compressed \code{cal.cbs} matrix. It makes quicker work
#' for Pareto/NBD estimation as implemented in this package, but the opposite is
#' true for BG/NBD. For proof, compare the definition of the
#' \code{\link{bgnbd.cbs.LL}} to that of \code{\link{pnbd.cbs.LL}}.
#'
#' @param params BG/NBD parameters - a vector with r, alpha, a, and b, in that
#' order. r and alpha are unobserved parameters for the NBD transaction
#' process. a and b are unobserved parameters for the Beta geometric dropout
#' process.
#' @param cal.cbs calibration period CBS (customer by sufficient statistic). It
#' must contain columns for frequency ("x"), recency ("t.x"), and total time
#' observed ("T.cal"). Note that recency must be the time between the start of
#' the calibration period and the customer's last transaction, not the time
#' between the customer's last transaction and the end of the calibration
#' period. If your data is compressed (see \code{\link{dc.compress.cbs}}),
#' a fourth column labeled "custs" (number of customers with a specific
#' combination of recency, frequency and length of calibration period) is
#' available.
#'
#' @seealso \code{\link{bgnbd.EstimateParameters}}
#' @seealso \code{\link{bgnbd.LL}}
#'
#' @return The total log-likelihood of the provided data.
#'
#' @examples
#' data(cdnowSummary)
#'
#' cal.cbs <- cdnowSummary$cbs
#' # cal.cbs already has column names required by method
#'
#' # random assignment of parameters
#' params <- c(0.5, 6, 1.2, 3.3)
#' # returns the log-likelihood of the given parameters
#' bgnbd.cbs.LL(params, cal.cbs)
#'
#' # compare the speed and results to the following:
#' cal.cbs.compressed <- dc.compress.cbs(cal.cbs)
#' bgnbd.cbs.LL (params, cal.cbs.compressed)
#'
#' # Returns the log likelihood of the parameters for a customer who
#' # made 3 transactions in a calibration period that ended at t=6,
#' # with the last transaction occurring at t=4.
#' bgnbd.LL(params, x=3, t.x=4, T.cal=6)
#'
#' # We can also give vectors as function parameters:
#' set.seed(7)
#' x <- sample(1:4, 10, replace = TRUE)
#' t.x <- sample(1:4, 10, replace = TRUE)
#' T.cal <- rep(4, 10)
#' bgnbd.LL(params, x, t.x, T.cal)
bgnbd.cbs.LL <- function(params,
cal.cbs) {
dc.check.model.params(printnames = c("r", "alpha", "a", "b"),
params = params,
func = "bgnbd.cbs.LL")
# Check that you have the right columns.
# They should be 'x', 't.x', 'T.cal' and optionally 'custs.'
# They stand for, respectively:
# -- x: frequency
# -- t.x: recency
# -- T.cal: observed calendar time
# -- custs: number of customers with this (x, t.x, T.cal) combo
foo <- colnames(cal.cbs)
stopifnot(all(c('x', 't.x', 'T.cal') %in% foo))
x <- cal.cbs[,'x']
t.x <- cal.cbs[,'t.x']
T.cal <- cal.cbs[,'T.cal']
# Avoid this unfurling exercise by calling bgnbd.cbs.LL
# with the uncompressed version of cal.cbs, which doesn't
# have a "custs" column.
if ("custs" %in% colnames(cal.cbs)) {
many_rows = function(vec, nreps) {
return(rep(1, nreps) %*% t.default(vec))
}
custs <- cal.cbs[, "custs"]
logvec = (1:length(custs)) * (custs > 1)
logvec = logvec[logvec > 0]
M = sum(logvec > 0)
for (i in 1:M) {
cal.cbs = rbind(cal.cbs,
many_rows(cal.cbs[logvec[i], ],
custs[logvec[i]] - 1))
}
x = cal.cbs[, "x"]
t.x = cal.cbs[, "t.x"]
T.cal = cal.cbs[, "T.cal"]
}
return(sum(bgnbd.LL(params, x, t.x, T.cal)))
}
#' BG/NBD Parameter Estimation
#'
#' Estimates parameters for the BG/NBD model.
#'
#' The best-fitting parameters are determined using the
#' \code{\link{bgnbd.cbs.LL}} function. The sum of the log-likelihood for each
#' customer (for a set of parameters) is maximized in order to estimate
#' parameters.
#'
#' A set of starting parameters must be provided for this method. If no
#' parameters are provided, (1,3,1,3) is used as a default. These values are
#' used because they provide good convergence across data sets. It may be useful
#' to use starting values for r and alpha that represent your best guess of the
#' heterogeneity in the buy and die rate of customers. It may be necessary to
#' run the estimation from multiple starting points to ensure that it converges.
#' To compare the log-likelihoods of different parameters, use
#' \code{\link{bgnbd.cbs.LL}}.
#'
#' The lower bound on the parameters to be estimated is always zero, since
#' BG/NBD parameters cannot be negative. The upper bound can be set with the
#' max.param.value parameter.
#'
#' This function may take some time to run.
#'
#' @param cal.cbs calibration period CBS (customer by sufficient statistic). It
#' must contain columns for frequency ("x"), recency ("t.x"), and total time
#' observed ("T.cal"). Note that recency must be the time between the start of
#' the calibration period and the customer's last transaction, not the time
#' between the customer's last transaction and the end of the calibration
#' period.
#' @param par.start initial BG/NBD parameters - a vector with r, alpha, a, and
#' b, in that order. r and alpha are unobserved parameters for the NBD
#' transaction process. a and b are unobserved parameters for the Beta
#' geometric dropout process.
#' @param max.param.value the upper bound on parameters.
#' @param method the optimization method(s) passed along to
#' \code{\link[optimx]{optimx}}.
#' @param hessian set it to TRUE if you want the Hessian matrix, and then you
#' might as well have the complete \code{\link[optimx]{optimx}} object
#' returned.
#' @return Vector of estimated parameters.
#' @seealso \code{\link{bgnbd.cbs.LL}}
#' @references Fader, Peter S.; Hardie, and Bruce G.S.. "Overcoming the BG/NBD
#' Model's #NUM! Error Problem." December. 2013. Web.
#' \url{http://brucehardie.com/notes/027/bgnbd_num_error.pdf}
#'
#' @examples
#' data(cdnowSummary)
#'
#' cal.cbs <- cdnowSummary$cbs
#' # cal.cbs already has column names required by method
#'
#' # starting-point parameters
#' startingparams <- c(1.0, 3, 1.0, 3)
#'
#' # estimated parameters
#' est.params <- bgnbd.EstimateParameters(cal.cbs = cal.cbs,
#' par.start = startingparams)
#'
#' # complete object returned by \code{\link[optimx]{optimx}}
#' optimx.set <- bgnbd.EstimateParameters(cal.cbs = cal.cbs,
#' par.start = startingparams,
#' hessian = TRUE)
#'
#' # log-likelihood of estimated parameters
#' bgnbd.cbs.LL(est.params, cal.cbs)
bgnbd.EstimateParameters <- function(cal.cbs,
par.start = c(1, 3, 1, 3),
max.param.value = 10000,
method = 'L-BFGS-B',
hessian = FALSE) {
dc.check.model.params(printnames = c("r", "alpha", "a", "b"),
params = par.start,
func = "bgnbd.EstimateParameters")
bgnbd.eLL <- function(params, cal.cbs, max.param.value) {
params <- exp(params)
params[params > max.param.value] = max.param.value
return(-1 * bgnbd.cbs.LL(params, cal.cbs))
}
logparams = log(par.start)
results <- optimx(par = logparams,
fn = bgnbd.eLL,
cal.cbs = cal.cbs,
max.param.value = max.param.value,
method = method,
hessian = hessian)
if(hessian == TRUE) {
message('Your parameter estimates are now on a log scale. Exponentiate them before use.')
return(results)
}
unlist(exp(results[method, c('p1', 'p2', 'p3', 'p4')]))
}
#' BG/NBD Conditional Expected Transactions
#'
#' E\[X(T.cal, T.cal + T.star) | x, t.x, r, alpha, a, b\]
#'
#' \code{T.star}, \code{x}, \code{t.x} and \code{T.cal} may be vectors. The
#' standard rules for vector operations apply - if they are not of the same
#' length, shorter vectors will be recycled (start over at the first element)
#' until they are as long as the longest vector. It is advisable to keep vectors
#' to the same length and to use single values for parameters that are to be the
#' same for all calculations. If one of these parameters has a length greater
#' than one, the output will be a vector of probabilities.
#'
#' @inheritParams bgnbd.LL
#' @param T.star length of time for which we are calculating the expected number
#' of transactions.
#' @param hardie if TRUE, use \code{\link{h2f1}} instead of
#' \code{\link[hypergeo]{hypergeo}}.
#' @return Number of transactions a customer is expected to make in a time
#' period of length t, conditional on their past behavior. If any of the input
#' parameters has a length greater than 1, this will be a vector of expected
#' number of transactions.
#' @seealso \code{\link{bgnbd.Expectation}}
#' @references Fader, Peter S.; Hardie, Bruce G.S.and Lee, Ka Lok. “Computing
#' P(alive) Using the BG/NBD Model.” December. 2008. Web.
#' \url{http://www.brucehardie.com/notes/021/palive_for_BGNBD.pdf}
#' @examples
#' params <- c(0.243, 4.414, 0.793, 2.426)
#' # Number of transactions a customer is expected to make in 2 time
#' # intervals, given that they made 10 repeat transactions in a time period
#' # of 39 intervals, with the 10th repeat transaction occurring in the 35th
#' # interval.
#' bgnbd.ConditionalExpectedTransactions(params, T.star=2, x=10, t.x=35, T.cal=39)
#'
#' # We can also compare expected transactions across different
#' # calibration period behaviors:
#' bgnbd.ConditionalExpectedTransactions(params, T.star=2, x=5:20, t.x=25, T.cal=39)
bgnbd.ConditionalExpectedTransactions <- function(params,
T.star,
x,
t.x,
T.cal,
hardie = TRUE) {
bgnbd.generalParams(params = params,
func = 'bgnbd.ConditionalExpectedTransactions',
x = x,
t.x = t.x,
T.cal = T.cal,
T.star = T.star,
hardie = hardie)$CET
}
#' BG/NBD Expectation
#'
#' Returns the number of repeat transactions that a randomly chosen customer
#' (for whom we have no prior information) is expected to make in a given time
#' period.
#'
#' E(X(t) | r, alpha, a, b)
#'
#' @param params BG/NBD parameters - a vector with r, alpha, a, and b, in that
#' order. r and alpha are unobserved parameters for the NBD transaction
#' process. a and b are unobserved parameters for the Beta geometric dropout
#' process.
#' @param t length of time for which we are calculating the expected number of
#' repeat transactions.
#' @param hardie if TRUE, use \code{\link{h2f1}} instead of
#' \code{\link[hypergeo]{hypergeo}}.
#' @return Number of repeat transactions a customer is expected to make in a
#' time period of length t.
#' @seealso \code{\link{bgnbd.ConditionalExpectedTransactions}}
#' @references Fader, Peter S.; Hardie, Bruce G.S.and Lee, Ka Lok. “Computing
#' P(alive) Using the BG/NBD Model.” December. 2008. Web.
#' \url{http://www.brucehardie.com/notes/021/palive_for_BGNBD.pdf}
#' @examples
#' params <- c(0.243, 4.414, 0.793, 2.426)
#'
#' # Number of repeat transactions a customer is expected to make in 2 time intervals.
#' bgnbd.Expectation(params, t=2, hardie = FALSE)
#'
#' # We can also compare expected transactions over time:
#' bgnbd.Expectation(params, t=1:10)
bgnbd.Expectation <- function(params,
t,
hardie = TRUE) {
dc.check.model.params(printnames = c("r", "alpha", "a", "b"),
params = params,
func = "bgnbd.Expectation")
if (any(t < 0) || !is.numeric(t))
stop("t must be numeric and may not contain negative numbers.")
r = params[1]
alpha = params[2]
a = params[3]
b = params[4]
term1 = (a + b - 1)/(a - 1)
term2 = (alpha/(alpha + t))^r
if(hardie == TRUE) {
term3 = h2f1(r, b, a + b - 1, t/(alpha + t))
} else {
term3 = Re(hypergeo(r, b, a + b - 1, t/(alpha + t)))
}
output = term1 * (1 - term2 * term3)
return(output)
}
#' BG/NBD Probability Mass Function
#'
#' Probability mass function for the BG/NBD.
#'
#' P(X(t)=x | r, alpha, a, b). Returns the probability that a customer makes x
#' repeat transactions in the time interval (0, t].
#'
#' Parameters t and x may be vectors. The standard rules for vector operations
#' apply - if they are not of the same length, the shorter vector will be
#' recycled (start over at the first element) until it is as long as the longest
#' vector. It is advisable to keep vectors to the same length and to use single
#' values for parameters that are to be the same for all calculations. If one of
#' these parameters has a length greater than one, the output will be a vector
#' of probabilities.
#'
#' @param params BG/NBD parameters - a vector with r, alpha, a, and b, in that
#' order. r and alpha are unobserved parameters for the NBD transaction
#' process. a and b are unobserved parameters for the Beta geometric dropout
#' process.
#' @param t length end of time period for which probability is being computed.
#' May also be a vector.
#' @param x number of repeat transactions by a random customer in the period
#' defined by t. May also be a vector.
#' @return Probability of X(t)=x conditional on model parameters. If t and/or x
#' has a length greater than one, a vector of probabilities will be returned.
#' @references Fader, Peter S.; Hardie, Bruce G.S.and Lee, Ka Lok. “Computing
#' P(alive) Using the BG/NBD Model.” December. 2008.
#' [Web.](http://www.brucehardie.com/notes/021/palive_for_BGNBD.pdf)
#' @examples
#' params <- c(0.243, 4.414, 0.793, 2.426)
#' # probability that a customer will make 10 repeat transactions in the
#' # time interval (0,2]
#' bgnbd.pmf(params, t=2, x=10)
#' # probability that a customer will make no repeat transactions in the
#' # time interval (0,39]
#' bgnbd.pmf(params, t=39, x=0)
#'
#' # Vectors may also be used as arguments:
#' bgnbd.pmf(params, t=30, x=11:20)
#' @md
bgnbd.pmf <- function(params,
t,
x) {
inputs <- try(dc.InputCheck(params = params,
func = 'bgnbd.pmf',
printnames = c("r", "alpha", "a", "b"),
x = x,
t = t))
if('try-error' == class(inputs)) return(inputs)
return(bgnbd.pmf.General(params,
t.start = 0,
t.end = inputs$t,
x = inputs$x))
}
#' Generalized BG/NBD Probability Mass Function
#'
#' Generalized probability mass function for the BG/NBD.
#'
#' P(X(t.start, t.end)=x | r, alpha, a, b). Returns the probability that a
#' customer makes x repeat transactions in the time interval (t.start, t.end\].
#'
#' It is impossible for a customer to make a negative number of repeat
#' transactions. This function will return an error if it is given negative
#' times or a negative number of repeat transactions. This function will also
#' return an error if t.end is less than t.start.
#'
#' t.start, t.end, and x may be vectors. The standard rules for vector
#' operations apply - if they are not of the same length, shorter vectors will
#' be recycled (start over at the first element) until they are as long as the
#' longest vector. It is advisable to keep vectors to the same length and to use
#' single values for parameters that are to be the same for all calculations. If
#' one of these parameters has a length greater than one, the output will be a
#' vector of probabilities.
#'
#' @param params BG/NBD parameters - a vector with r, alpha, a, and b, in that
#' order. r and alpha are unobserved parameters for the NBD transaction
#' process. a and b are unobserved parameters for the Beta geometric dropout
#' process.
#' @param t.start start of time period for which probability is being
#' calculated. It can also be a vector of values.
#' @param t.end end of time period for which probability is being calculated.
#' It can also be a vector of values.
#' @param x number of repeat transactions by a random customer in the period
#' defined by (t.start, t.end]. It can also be a vector of values.
#' @return Probability of x transaction occuring between t.start and t.end
#' conditional on model parameters. If t.start, t.end, and/or x has a length
#' greater than one, a vector of probabilities will be returned.
#' @references Fader, Peter S.; Hardie, Bruce G.S.and Lee, Ka Lok. “Computing
#' P(alive) Using the BG/NBD Model.” December. 2008.
#' [Web.](http://www.brucehardie.com/notes/021/palive_for_BGNBD.pdf)
#' @examples
#' params <- c(0.243, 4.414, 0.793, 2.426)
#' # probability that a customer will make 10 repeat transactions in the
#' # time interval (1,2]
#' bgnbd.pmf.General(params, t.start=1, t.end=2, x=10)
#' # probability that a customer will make no repeat transactions in the
#' # time interval (39,78]
#' bgnbd.pmf.General(params, t.start=39, t.end=78, x=0)
#' @md
bgnbd.pmf.General <- function(params,
t.start,
t.end,
x) {
inputs <- try(dc.InputCheck(params = params,
func = 'bgnbd.pmf.General',
printnames = c("r", "alpha", "a", "b"),
t.start = t.start,
t.end = t.end,
x = x))
if('try-error' == class(inputs)) return(inputs)
t.start = inputs$t.start
t.end = inputs$t.end
x = inputs$x
max.length <- nrow(inputs)
if (any(t.start > t.end)) {
stop("Error in bgnbd.pmf.General: t.start > t.end.")
}
r <- params[1]
alpha <- params[2]
a <- params[3]
b <- params[4]
equation.part.0 <- rep(0, max.length)
t = t.end - t.start
term3 = rep(0, max.length)
term1 = beta(a, b + x)/beta(a, b) *
gamma(r + x)/gamma(r)/factorial(x) *
((alpha/(alpha + t))^r) * ((t/(alpha + t))^x)
for (i in 1:max.length) {
if (x[i] > 0) {
ii = c(0:(x[i] - 1))
summation.term = sum(gamma(r + ii)/gamma(r)/factorial(ii) *
((t[i]/(alpha + t[i]))^ii))
term3[i] = 1 - (((alpha/(alpha + t[i]))^r) * summation.term)
}
}
term2 = as.numeric(x > 0) * beta(a + 1, b + x - 1)/beta(a, b) * term3
return(term1 + term2)
}
#' BG/NBD P(Alive)
#'
#' Uses BG/NBD model parameters and a customer's past transaction behavior to
#' return the probability that they are still alive at the end of the
#' calibration period.
#'
#' P(Alive | X=x, t.x, T.cal, r, alpha, a, b)
#'
#' x, t.x, and T.cal may be vectors. The standard rules for vector operations
#' apply - if they are not of the same length, shorter vectors will be recycled
#' (start over at the first element) until they are as long as the longest
#' vector. It is advisable to keep vectors to the same length and to use single
#' values for parameters that are to be the same for all calculations. If one of
#' these parameters has a length greater than one, the output will be a vector
#' of probabilities.
#'
#' @inheritParams bgnbd.LL
#' @return Probability that the customer is still alive at the end of the
#' calibration period. If x, t.x, and/or T.cal has a length greater than one,
#' then this will be a vector of probabilities (containing one element
#' matching each element of the longest input vector).
#' @references Fader, Peter S.; Hardie, Bruce G.S.and Lee, Ka Lok. “Computing
#' P(alive) Using the BG/NBD Model.” December. 2008.
#' [Web.](http://www.brucehardie.com/notes/021/palive_for_BGNBD.pdf)
#' @examples
#' params <- c(0.243, 4.414, 0.793, 2.426)
#'
#' bgnbd.PAlive(params, x=23, t.x=39, T.cal=39)
#' # P(Alive) of a customer who has the same recency and total
#' # time observed.
#'
#' bgnbd.PAlive(params, x=5:20, t.x=30, T.cal=39)
#' # Note the "increasing frequency paradox".
#'
#' # To visualize the distribution of P(Alive) across customers:
#'
#' data(cdnowSummary)
#' cbs <- cdnowSummary$cbs
#' params <- bgnbd.EstimateParameters(cbs, par.start = c(0.243, 4.414, 0.793, 2.426))
#' p.alives <- bgnbd.PAlive(params, cbs[,"x"], cbs[,"t.x"], cbs[,"T.cal"])
#' plot(density(p.alives))
#' @md
bgnbd.PAlive <- function(params,
x,
t.x,
T.cal) {
bgnbd.generalParams(params = params,
func = 'bgnbd.PAlive',
x = x,
t.x = t.x,
T.cal = T.cal)$PAlive
}
#' BG/NBD Expected Cumulative Transactions
#'
#' Calculates the expected cumulative total repeat transactions by all customers
#' for the calibration and holdout periods.
#'
#' The function automatically divides the total period up into n.periods.final
#' time intervals. n.periods.final does not have to be in the same unit of time
#' as the T.cal data. For example: - if your T.cal data is in weeks, and you
#' want cumulative transactions per week, n.periods.final would equal T.star. -
#' if your T.cal data is in weeks, and you want cumulative transactions per day,
#' n.periods.final would equal T.star * 7.
#'
#' The holdout period should immediately follow the calibration period. This
#' function assume that all customers' calibration periods end on the same date,
#' rather than starting on the same date (thus customers' birth periods are
#' determined using max(T.cal) - T.cal rather than assuming that it is 0).
#'
#' @param params BG/NBD parameters - a vector with r, alpha, a, and b, in that
#' order. r and alpha are unobserved parameters for the NBD transaction
#' process. a and b are unobserved parameters for the Beta geometric dropout
#' process.
#' @param T.cal a vector to represent customers' calibration period lengths
#' (in other words, the "T.cal" column from a customer-by-sufficient-statistic
#' matrix).
#' @param T.tot end of holdout period. Must be a single value, not a vector.
#' @param n.periods.final number of time periods in the calibration and
#' holdout periods. See details.
#' @param hardie if TRUE, use h2f1 instead of hypergeo.
#' @return Vector of expected cumulative total repeat transactions by all
#' customers.
#' @seealso [`bgnbd.Expectation`]
#' @examples
#' data(cdnowSummary)
#'
#' cal.cbs <- cdnowSummary$cbs
#' # cal.cbs already has column names required by method
#'
#' params <- c(0.243, 4.414, 0.793, 2.426)
#'
#' # Returns a vector containing cumulative repeat transactions for 273 days.
#' # All parameters are in weeks; the calibration period lasted 39 weeks.
#' bgnbd.ExpectedCumulativeTransactions(params,
#' T.cal = cal.cbs[,"T.cal"],
#' T.tot = 39,
#' n.periods.final = 273,
#' hardie = TRUE)
#' @md
bgnbd.ExpectedCumulativeTransactions <- function(params,
T.cal,
T.tot,
n.periods.final,
hardie = TRUE) {
dc.check.model.params(printnames = c("r", "alpha", "s", "beta"),
params = params,
func = "bgnbd.ExpectedCumulativeTransactions")
if (any(T.cal < 0) || !is.numeric(T.cal))
stop("T.cal must be numeric and may not contain negative numbers.")
if (length(T.tot) > 1 || T.tot < 0 || !is.numeric(T.tot))
stop("T.cal must be a single numeric value and may not be negative.")
if (length(n.periods.final) > 1 || n.periods.final < 0 || !is.numeric(n.periods.final))
stop("n.periods.final must be a single numeric value and may not be negative.")
intervals <- seq(T.tot/n.periods.final,
T.tot,
length.out = n.periods.final)
cust.birth.periods <- max(T.cal) - T.cal
expected.transactions <- sapply(intervals,
function(interval) {
if (interval <= min(cust.birth.periods)) return(0)
t <- interval - cust.birth.periods[cust.birth.periods <= interval]
sum(bgnbd.Expectation(params = params,
t = t,
hardie = hardie))
})
return(expected.transactions)
}
#' BG/NBD Plot Frequency in Calibration Period
#'
#' Plots a histogram and returns a matrix comparing the actual and expected
#' number of customers who made a certain number of repeat transactions in the
#' calibration period, binned according to calibration period frequencies.
#'
#' This function requires a censor number, which cannot be higher than the
#' highest frequency in the calibration period CBS. The output matrix will have
#' (censor + 1) bins, starting at frequencies of 0 transactions and ending at a
#' bin representing calibration period frequencies at or greater than the censor
#' number. The plot may or may not include a bin for zero frequencies, depending
#' on the plotZero parameter.
#'
#' @param params BG/NBD parameters - a vector with r, alpha, a, and b, in that
#' order. r and alpha are unobserved parameters for the NBD transaction
#' process. a and b are unobserved parameters for the Beta geometric dropout
#' process.
#' @param cal.cbs calibration period CBS (customer by sufficient statistic). It
#' must contain columns for frequency ("x") and total time observed ("T.cal").
#' @param censor integer used to censor the data. See details.
#' @param plotZero If FALSE, the histogram will exclude the zero bin.
#' @param xlab descriptive label for the x axis.
#' @param ylab descriptive label for the y axis.
#' @param title title placed on the top-center of the plot.
#' @return Calibration period repeat transaction frequency comparison matrix
#' (actual vs. expected).
#' @examples
#' data(cdnowSummary)
#'
#' cal.cbs <- cdnowSummary$cbs
#' # cal.cbs already has column names required by method
#'
#' # parameters estimated using bgnbd.EstimateParameters
#' est.params <- c(0.243, 4.414, 0.793, 2.426)
#' # the maximum censor number that can be used
#' max(cal.cbs[,"x"])
#'
#' bgnbd.PlotFrequencyInCalibration(est.params, cal.cbs, censor=7)
bgnbd.PlotFrequencyInCalibration <- function(params,
cal.cbs,
censor,
plotZero = TRUE,
xlab = "Calibration period transactions",
ylab = "Customers",
title = "Frequency of Repeat Transactions") {
tryCatch(x <- cal.cbs[, "x"], error = function(e) stop("Error in bgnbd.PlotFrequencyInCalibration: cal.cbs must have a frequency column labelled \"x\""))
tryCatch(T.cal <- cal.cbs[, "T.cal"], error = function(e) stop("Error in bgnbd.PlotFrequencyInCalibration: cal.cbs must have a column for length of time observed labelled \"T.cal\""))
dc.check.model.params(c("r", "alpha", "a", "b"), params, "bgnbd.PlotFrequencyInCalibration")
if (censor > max(x))
stop("censor too big (> max freq) in PlotFrequencyInCalibration.")
x = cal.cbs[, "x"]
T.cal = cal.cbs[, "T.cal"]
n.x <- rep(0, max(x) + 1)
ncusts = nrow(cal.cbs)
for (ii in unique(x)) {
# Get number of customers to buy n.x times, over the grid of all possible n.x
# values (no censoring)
n.x[ii + 1] <- sum(ii == x)
}
n.x.censor <- sum(n.x[(censor + 1):length(n.x)])
n.x.actual <- c(n.x[1:censor], n.x.censor) # This upper truncates at censor (ie. if censor=7, 8 categories: {0, 1, ..., 6, 7+}).
T.value.counts <- table(T.cal) # This is the table of counts of all time durations from customer birth to end of calibration period.
T.values <- as.numeric(names(T.value.counts)) # These are all the unique time durations from customer birth to end of calibration period.
n.T.values <- length(T.values) # These are the number of time durations we need to consider.
n.x.expected <- rep(0, length(n.x.actual)) # We'll store the probabilities in here.
n.x.expected.all <- rep(0, max(x) + 1) # We'll store the probabilities in here.
for (ii in 0:max(x)) {
# We want to run over the probability of each transaction amount.
this.x.expected = 0
for (T.idx in 1:n.T.values) {
# We run over all people who had all time durations.
T = T.values[T.idx]
if (T == 0)
next
n.T = T.value.counts[T.idx] # This is the number of customers who had this time duration.
prob.of.this.x.for.this.T = bgnbd.pmf(params, T, ii)
expected.given.x.and.T = n.T * prob.of.this.x.for.this.T
this.x.expected = this.x.expected + expected.given.x.and.T
}
n.x.expected.all[ii + 1] = this.x.expected
}
n.x.expected[1:censor] = n.x.expected.all[1:censor]
n.x.expected[censor + 1] = sum(n.x.expected.all[(censor + 1):(max(x) + 1)])
col.names <- paste(rep("freq", length(censor + 1)), (0:censor), sep = ".")
col.names[censor + 1] <- paste(col.names[censor + 1], "+", sep = "")
censored.freq.comparison <- rbind(n.x.actual, n.x.expected)
colnames(censored.freq.comparison) <- col.names
cfc.plot <- censored.freq.comparison
if (plotZero == FALSE)
cfc.plot <- cfc.plot[, -1]
n.ticks <- ncol(cfc.plot)
if (plotZero == TRUE) {
x.labels <- 0:(n.ticks - 1)
x.labels[n.ticks] <- paste(n.ticks - 1, "+", sep = "")
}
ylim <- c(0, ceiling(max(cfc.plot) * 1.1))
barplot(cfc.plot, names.arg = x.labels, beside = TRUE, ylim = ylim, main = title,
xlab = xlab, ylab = ylab, col = 1:2)
legend("topright", legend = c("Actual", "Model"), col = 1:2, lwd = 2)
return(censored.freq.comparison)
}
#' BG/NBD Plot Frequency vs. Conditional Expected Frequency
#'
#' Plots the actual and conditional expected number transactions made by
#' customers in the holdout period, binned according to calibration period
#' frequencies. Also returns a matrix with this comparison and the number of
#' customers in each bin.
#'
#' This function requires a censor number, which cannot be higher than the
#' highest frequency in the calibration period CBS. The output matrix will have
#' (censor + 1) bins, starting at frequencies of 0 transactions and ending at a
#' bin representing calibration period frequencies at or greater than the censor
#' number.
#'
#' @param params BG/NBD parameters - a vector with r, alpha, a, and b, in that
#' order. r and alpha are unobserved parameters for the NBD transaction
#' process. a and b are unobserved parameters for the Beta geometric dropout
#' process.
#' @param T.star length of then holdout period.
#' @param cal.cbs calibration period CBS (customer by sufficient statistic).
#' It must contain columns for frequency ("x"), recency ("t.x"), and total
#' time observed ("T.cal"). Note that recency must be the time between the
#' start of the calibration period and the customer's last transaction, not
#' the time between the customer's last transaction and the end of the
#' calibration period.
#' @param x.star vector of transactions made by each customer in the holdout
#' period.
#' @param censor integer used to censor the data. See details.
#' @param xlab descriptive label for the x axis.
#' @param ylab descriptive label for the y axis.
#' @param xticklab vector containing a label for each tick mark on the x axis.
#' @param title title placed on the top-center of the plot.
#' @return Holdout period transaction frequency comparison matrix (actual vs.
#' expected).
#' @examples
#' data(cdnowSummary)
#'
#' cal.cbs <- cdnowSummary$cbs
#' # cal.cbs already has column names required by method
#'
#' # number of transactions by each customer in the 39 weeks
#' # following the calibration period
#' x.star <- cal.cbs[,"x.star"]
#'
#' # parameters estimated using bgnbd.EstimateParameters
#' est.params <- c(0.243, 4.414, 0.793, 2.426)
#' # the maximum censor number that can be used
#' max(cal.cbs[,"x"])
#'
#' # plot conditional expected holdout period frequencies,
#' # binned according to calibration period frequencies
#' bgnbd.PlotFreqVsConditionalExpectedFrequency(est.params,
#' T.star = 39,
#' cal.cbs,
#' x.star,
#' censor = 7)
bgnbd.PlotFreqVsConditionalExpectedFrequency <- function(params,
T.star,
cal.cbs,
x.star,
censor,
xlab = "Calibration period transactions",
ylab = "Holdout period transactions",
xticklab = NULL,
title = "Conditional Expectation") {
tryCatch(x <- cal.cbs[, "x"],
error = function(e) stop("Error in bgnbd.PlotFreqVsConditionalExpectedFrequency: cal.cbs must have a frequency column labelled \"x\""))
tryCatch(t.x <- cal.cbs[, "t.x"],
error = function(e) stop("Error in bgnbd.PlotFreqVsConditionalExpectedFrequency: cal.cbs must have a recency column labelled \"t.x\""))
tryCatch(T.cal <- cal.cbs[, "T.cal"],
error = function(e) stop("Error in bgnbd.PlotFreqVsConditionalExpectedFrequency: cal.cbs must have a column for length of time observed labelled \"T.cal\""))
dc.check.model.params(c("r", "alpha", "a", "b"), params, "bgnbd.PlotFreqVsConditionalExpectedFrequency")
if (censor > max(x))
stop("censor too big (> max freq) in PlotFreqVsConditionalExpectedFrequency.")
if (any(T.star < 0) || !is.numeric(T.star))
stop("T.star must be numeric and may not contain negative numbers.")
if (any(x.star < 0) || !is.numeric(x.star))
stop("x.star must be numeric and may not contain negative numbers.")
n.bins = censor + 1
transaction.actual = rep(0, n.bins)
transaction.expected = rep(0, n.bins)
bin.size = rep(0, n.bins)
for (cc in 0:censor) {
if (cc != censor) {
this.bin = which(cc == x)
} else if (cc == censor) {
this.bin = which(x >= cc)
}
n.this.bin = length(this.bin)
bin.size[cc + 1] = n.this.bin
transaction.actual[cc + 1] = sum(x.star[this.bin])/n.this.bin
transaction.expected[cc + 1] = sum(bgnbd.ConditionalExpectedTransactions(params,
T.star, x[this.bin], t.x[this.bin], T.cal[this.bin]))/n.this.bin
}
col.names = paste(rep("freq", length(censor + 1)), (0:censor), sep = ".")
col.names[censor + 1] = paste(col.names[censor + 1], "+", sep = "")
comparison = rbind(transaction.actual, transaction.expected, bin.size)
colnames(comparison) = col.names
if (is.null(xticklab) == FALSE) {
x.labels = xticklab
}
if (is.null(xticklab) != FALSE) {
if (censor < ncol(comparison)) {
x.labels = 0:(censor)
x.labels[censor + 1] = paste(censor, "+", sep = "")
}
if (censor >= ncol(comparison)) {
x.labels = 0:(ncol(comparison))
}
}
actual = comparison[1, ]
expected = comparison[2, ]
ylim = c(0, ceiling(max(c(actual, expected)) * 1.1))
plot(actual, type = "l", xaxt = "n", col = 1, ylim = ylim, xlab = xlab, ylab = ylab,
main = title)
lines(expected, lty = 2, col = 2)
axis(1, at = 1:ncol(comparison), labels = x.labels)
legend("topleft", legend = c("Actual", "Model"), col = 1:2, lty = 1:2, lwd = 1)
return(comparison)
}
#' BG/NBD Plot Actual vs. Conditional Expected Frequency by Recency
#'
#' Plots the actual and conditional expected number of transactions made by
#' customers in the holdout period, binned according to calibration period
#' recencies. Also returns a matrix with this comparison and the number of
#' customers in each bin.
#'
#' This function does bin customers exactly according to recency; it bins
#' customers according to integer units of the time period of cal.cbs.
#' Therefore, if you are using weeks in your data, customers will be binned as
#' follows: customers with recencies between the start of the calibration period
#' (inclusive) and the end of week one (exclusive); customers with recencies
#' between the end of week one (inclusive) and the end of week two (exclusive);
#' etc.
#'
#' The matrix and plot will contain the actual number of transactions made by
#' each bin in the holdout period, as well as the expected number of
#' transactions made by that bin in the holdout period, conditional on that
#' bin's behavior during the calibration period.
#'
#' @inheritParams bgnbd.PlotFreqVsConditionalExpectedFrequency
#' @return Matrix comparing actual and conditional expected transactions in the
#' holdout period.
#' @examples
#' data(cdnowSummary)
#'
#' cal.cbs <- cdnowSummary$cbs
#' # cal.cbs already has column names required by method
#'
#' # number of transactions by each customer in the 39 weeks following
#' # the calibration period
#' x.star <- cal.cbs[,"x.star"]
#'
#' # parameters estimated using bgnbd.EstimateParameters
#' est.params <- c(0.243, 4.414, 0.793, 2.426)
#'
#' # plot conditional expected holdout period transactions,
#' # binned according to calibration period recencies
#' bgnbd.PlotRecVsConditionalExpectedFrequency(est.params,
#' cal.cbs,
#' T.star = 39,
#' x.star)
bgnbd.PlotRecVsConditionalExpectedFrequency <- function(params,
cal.cbs,
T.star,
x.star,
xlab = "Calibration period recency",
ylab = "Holdout period transactions",
xticklab = NULL,
title = "Actual vs. Conditional Expected Transactions by Recency") {
dc.check.model.params(c("r", "alpha", "a", "b"), params, "bgnbd.PlotRecVsConditionalExpectedFrequency")
if (any(T.star < 0) || !is.numeric(T.star))
stop("T.star must be numeric and may not contain negative numbers.")
if (any(x.star < 0) || !is.numeric(x.star))
stop("x.star must be numeric and may not contain negative numbers.")
tryCatch(x <- cal.cbs[, "x"],
error = function(e) stop("Error in bgnbd.PlotRecVsConditionalExpectedFrequency: cal.cbs must have a frequency column labelled \"x\""))
tryCatch(t.x <- cal.cbs[, "t.x"],
error = function(e) stop("Error in bgnbd.PlotRecVsConditionalExpectedFrequency: cal.cbs must have a recency column labelled \"t.x\""))
tryCatch(T.cal <- cal.cbs[, "T.cal"],
error = function(e) stop("Error in bgnbd.PlotRecVsConditionalExpectedFrequency: cal.cbs must have a column for length of time observed labelled \"T.cal\""))
t.values <- sort(unique(t.x))
n.recs <- length(t.values)
transaction.actual <- rep(0, n.recs)
transaction.expected <- rep(0, n.recs)
rec.size <- rep(0, n.recs)
for (tt in 1:n.recs) {
this.t.x <- t.values[tt]
this.rec <- which(t.x == this.t.x)
n.this.rec <- length(this.rec)
rec.size[tt] <- n.this.rec
transaction.actual[tt] <- sum(x.star[this.rec])/n.this.rec
transaction.expected[tt] <- sum(bgnbd.ConditionalExpectedTransactions(params,
T.star, x[this.rec], t.x[this.rec], T.cal[this.rec]))/n.this.rec
}
comparison <- rbind(transaction.actual, transaction.expected, rec.size)
colnames(comparison) <- round(t.values, 3)
bins <- seq(1, ceiling(max(t.x)))
n.bins <- length(bins)
actual <- rep(0, n.bins)
expected <- rep(0, n.bins)
bin.size <- rep(0, n.bins)
x.labels <- NULL
if (is.null(xticklab) == FALSE) {
x.labels <- xticklab
} else {
x.labels <- 1:(n.bins)
}
point.labels <- rep("", n.bins)
point.y.val <- rep(0, n.bins)
for (ii in 1:n.bins) {
if (ii < n.bins) {
this.bin <- which(as.numeric(colnames(comparison)) >= (ii - 1) & as.numeric(colnames(comparison)) <
ii)
} else if (ii == n.bins) {
this.bin <- which(as.numeric(colnames(comparison)) >= ii - 1)
}
actual[ii] <- sum(comparison[1, this.bin])/length(comparison[1, this.bin])
expected[ii] <- sum(comparison[2, this.bin])/length(comparison[2, this.bin])
bin.size[ii] <- sum(comparison[3, this.bin])
}
ylim <- c(0, ceiling(max(c(actual, expected)) * 1.1))
plot(actual, type = "l", xaxt = "n", col = 1, ylim = ylim, xlab = xlab, ylab = ylab,
main = title)
lines(expected, lty = 2, col = 2)
axis(1, at = 1:n.bins, labels = x.labels)
legend("topleft", legend = c("Actual", "Model"), col = 1:2, lty = 1:2, lwd = 1)
return(rbind(actual, expected, bin.size))
}
#' BG/NBD Plot Transaction Rate Heterogeneity
#'
#' Plots and returns the estimated gamma distribution of lambda (customers'
#' propensities to purchase).
#'
#' This returns the distribution of each customer's Poisson parameter, which
#' determines the rate at which each customer buys.
#'
#' @param params BG/NBD parameters - a vector with r, alpha, a, and b, in that
#' order. r and alpha are unobserved parameters for the NBD transaction
#' process. a and b are unobserved parameters for the Beta geometric dropout
#' process.
#' @param lim upper-bound of the x-axis. A number is chosen by the function if
#' none is provided.
#' @return Distribution of customers' propensities to purchase.
#' @examples
#' params <- c(0.243, 4.414, 0.793, 2.426)
#' bgnbd.PlotTransactionRateHeterogeneity(params)
#' params <- c(0.53, 4.414, 0.793, 2.426)
#' bgnbd.PlotTransactionRateHeterogeneity(params)
bgnbd.PlotTransactionRateHeterogeneity <- function(params,
lim = NULL) {
dc.check.model.params(c("r", "alpha", "a", "b"), params, "bgnbd.PlotTransactionRateHeterogeneity")
shape <- params[1]
rate <- params[2]
rate.mean <- round(shape/rate, 4)
rate.var <- round(shape/rate^2, 4)
if (is.null(lim)) {
lim = qgamma(0.99, shape = shape, rate = rate)
}
x.axis.ticks <- seq(0, lim, length.out = 100)
heterogeneity <- dgamma(x.axis.ticks, shape = shape, rate = rate)
plot(x.axis.ticks, heterogeneity, type = "l", xlab = "Transaction Rate", ylab = "Density",
main = "Heterogeneity in Transaction Rate")
mean.var.label <- paste("Mean:", rate.mean, " Var:", rate.var)
mtext(mean.var.label, side = 3)
return(rbind(x.axis.ticks, heterogeneity))
}
#' BG/NBD Plot Dropout Probability Heterogeneity
#'
#' Plots and returns the estimated gamma distribution of p (customers'
#' probability of dropping out immediately after a transaction).
#'
#' @inheritParams bgnbd.PlotTransactionRateHeterogeneity
#' @return Distribution of customers' probabilities of dropping out.
#' @examples
#' params <- c(0.243, 4.414, 0.793, 2.426)
#' bgnbd.PlotDropoutRateHeterogeneity(params)
#' params <- c(0.243, 4.414, 1.33, 2.426)
#' bgnbd.PlotDropoutRateHeterogeneity(params)
bgnbd.PlotDropoutRateHeterogeneity <- function(params,
lim = NULL) {
dc.check.model.params(c("r", "alpha", "a", "b"), params, "bgnbd.PlotDropoutRateHeterogeneity")
alpha_param = params[3]
beta_param = params[4]
beta_param.mean = round(alpha_param/(alpha_param + beta_param), 4)
beta_param.var = round(alpha_param * beta_param/((alpha_param + beta_param)^2)/(alpha_param +
beta_param + 1), 4)
if (is.null(lim)) {
# get right end point of grid
lim = qbeta(0.99, shape1 = alpha_param, shape2 = beta_param)
}
x.axis.ticks = seq(0, lim, length.out = 100)
heterogeneity = dbeta(x.axis.ticks, shape1 = alpha_param, shape2 = beta_param)
plot(x.axis.ticks, heterogeneity, type = "l", xlab = "Dropout Probability p",
ylab = "Density", main = "Heterogeneity in Dropout Probability")
mean.var.label = paste("Mean:", beta_param.mean, " Var:", beta_param.var)
mtext(mean.var.label, side = 3)
return(rbind(x.axis.ticks, heterogeneity))
}
#' BG/NBD Tracking Cumulative Transactions Plot
#'
#' Plots the actual and expected cumulative total repeat transactions by all
#' customers for the calibration and holdout periods, and returns this
#' comparison in a matrix.
#'
#' actual.cu.tracking.data does not have to be in the same unit of time as the
#' T.cal data. T.tot will automatically be divided into periods to match the
#' length of actual.cu.tracking.data. See
#' [bgnbd.ExpectedCumulativeTransactions].
#'
#' The holdout period should immediately follow the calibration period. This
#' function assume that all customers' calibration periods end on the same date,
#' rather than starting on the same date (thus customers' birth periods are
#' determined using max(T.cal) - T.cal rather than assuming that it is 0).
#'
#' @inheritParams bgnbd.ExpectedCumulativeTransactions
#' @inheritParams bgnbd.PlotFreqVsConditionalExpectedFrequency
#' @param actual.cu.tracking.data vector containing the cumulative number of
#' repeat transactions made by customers for each period in the total time
#' period (both calibration and holdout periods). See details.
#' @return Matrix containing actual and expected cumulative repeat transactions.
#' @examples
#' data(cdnowSummary)
#'
#' cal.cbs <- cdnowSummary$cbs
#' # cal.cbs already has column names required by method
#'
#' # Cumulative repeat transactions made by all customers across calibration
#' # and holdout periods
#' cu.tracking <- cdnowSummary$cu.tracking
#'
#' # parameters estimated using bgnbd.EstimateParameters
#' est.params <- c(0.243, 4.414, 0.793, 2.426)
#'
#' # All parameters are in weeks; the calibration period lasted 39
#' # weeks and the holdout period another 39.
#' bgnbd.PlotTrackingCum(est.params,
#' T.cal = cal.cbs[,"T.cal"],
#' T.tot = 78,
#' actual.cu.tracking.data = cu.tracking,
#' hardie = TRUE)
#' @md
bgnbd.PlotTrackingCum <- function(params,
T.cal,
T.tot,
actual.cu.tracking.data,
n.periods.final = NA,
hardie = TRUE,
xlab = "Week",
ylab = "Cumulative Transactions",
xticklab = NULL,
title = "Tracking Cumulative Transactions") {
dc.check.model.params(c("r", "alpha", "a", "b"), params, "bgnbd.Plot.PlotTrackingCum")
if (any(T.cal < 0) || !is.numeric(T.cal))
stop("T.cal must be numeric and may not contain negative numbers.")
if (any(actual.cu.tracking.data < 0) || !is.numeric(actual.cu.tracking.data))
stop("actual.cu.tracking.data must be numeric and may not contain negative numbers.")
if (length(T.tot) > 1 || T.tot < 0 || !is.numeric(T.tot))
stop("T.cal must be a single numeric value and may not be negative.")
actual <- actual.cu.tracking.data
if(is.na(n.periods.final)) n.periods.final <- length(actual)
expected <- bgnbd.ExpectedCumulativeTransactions(params,
T.cal,
T.tot,
n.periods.final,
hardie)
cu.tracking.comparison <- rbind(actual, expected)
ylim <- c(0, max(c(actual, expected)) * 1.05)
plot(actual, type = "l", xaxt = "n", xlab = xlab, ylab = ylab, col = 1, ylim = ylim,
main = title)
lines(expected, lty = 2, col = 2)
if (is.null(xticklab) == FALSE) {
if (ncol(cu.tracking.comparison) != length(xticklab)) {
stop("Plot error, xticklab does not have the correct size")
}
axis(1, at = 1:ncol(cu.tracking.comparison), labels = xticklab)
} else {
axis(1, at = 1:length(actual), labels = 1:length(actual))
}
abline(v = max(T.cal), lty = 2)
legend("bottomright", legend = c("Actual", "Model"), col = 1:2, lty = 1:2, lwd = 1)
return(cu.tracking.comparison)
}
#' BG/NBD Tracking Incremental Transactions Comparison
#'
#' Plots the actual and expected incremental total repeat transactions by all
#' customers for the calibration and holdout periods, and returns this
#' comparison in a matrix.
#'
#' actual.inc.tracking.data does not have to be in the same unit of time as the
#' T.cal data. T.tot will automatically be divided into periods to match the
#' length of actual.inc.tracking.data. See
#' [bgnbd.ExpectedCumulativeTransactions].
#'
#' The holdout period should immediately follow the calibration period. This
#' function assume that all customers' calibration periods end on the same date,
#' rather than starting on the same date (thus customers' birth periods are
#' determined using max(T.cal) - T.cal rather than assuming that it is 0).
#'
#' @inheritParams bgnbd.PlotTrackingCum
#' @param actual.inc.tracking.data vector containing the incremental number of
#' repeat transactions made by customers for each period in the total time
#' period (both calibration and holdout periods). See details.
#' @return Matrix containing actual and expected incremental repeat
#' transactions.
#' @examples
#' data(cdnowSummary)
#' cal.cbs <- cdnowSummary$cbs
#' # cal.cbs already has column names required by method
#'
#' # Cumulative repeat transactions made by all customers across calibration
#' # and holdout periods
#' cu.tracking <- cdnowSummary$cu.tracking
#' # make the tracking data incremental
#' inc.tracking <- dc.CumulativeToIncremental(cu.tracking)
#'
#' # parameters estimated using bgnbd.EstimateParameters
#' est.params <- c(0.243, 4.414, 0.793, 2.426)
#'
#' # All parameters are in weeks; the calibration period lasted 39
#' # weeks and the holdout period another 39.
#' bgnbd.PlotTrackingInc(est.params, ,
#' T.cal = cal.cbs[,"T.cal"],
#' T.tot = 78,
#' actual.inc.tracking.data = inc.tracking,
#' hardie = TRUE)
#' @md
bgnbd.PlotTrackingInc <- function(params,
T.cal,
T.tot,
actual.inc.tracking.data,
n.periods.final = NA,
hardie = TRUE,
xlab = "Week",
ylab = "Transactions",
xticklab = NULL,
title = "Tracking Weekly Transactions") {
dc.check.model.params(printnames = c("r", "alpha", "a", "b"),
params = params,
func = "bgnbd.Plot.PlotTrackingCum")
if (any(T.cal < 0) || !is.numeric(T.cal))
stop("T.cal must be numeric and may not contain negative numbers.")
if (any(actual.inc.tracking.data < 0) || !is.numeric(actual.inc.tracking.data))
stop("actual.inc.tracking.data must be numeric and may not contain negative numbers.")
if (length(T.tot) > 1 || T.tot < 0 || !is.numeric(T.tot))
stop("T.cal must be a single numeric value and may not be negative.")
actual <- actual.inc.tracking.data
if(is.na(n.periods.final)) n.periods.final <- length(actual)
expected <- dc.CumulativeToIncremental(bgnbd.ExpectedCumulativeTransactions(params,
T.cal,
T.tot,
n.periods.final,
hardie))
ylim <- c(0, max(c(actual, expected)) * 1.05)
plot(actual, type = "l", xaxt = "n", xlab = xlab, ylab = ylab, col = 1, ylim = ylim,
main = title)
lines(expected, lty = 2, col = 2)
if (is.null(xticklab) == FALSE) {
if (length(actual) != length(xticklab)) {
stop("Plot error, xticklab does not have the correct size")
}
axis(1, at = 1:length(actual), labels = xticklab)
} else {
axis(1, at = 1:length(actual), labels = 1:length(actual))
}
abline(v = max(T.cal), lty = 2)
legend("topright", legend = c("Actual", "Model"), col = 1:2, lty = 1:2, lwd = 1)
return(rbind(actual, expected))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.