R/BG.R

In sriharitn/foretell: Projecting Customer Retention Based on Fader and Hardie Probability Models

#' Beta Geometric (BG) Model for Projecting Customer Retention.
#'
#' \code{BG} is a beta geometric model implemented based on \code{Fader and Hardie} probability based projection methedology. The survivor function for \code{BG} is \deqn{Beta(a,b+t)/Beta(a,b)}
#'
#' @param surv_value a numeric vector of historical customer retention percentage should start at 100 and non-starting values should be between 0 and less than 100
#' @param h forecasting horizon
#' @param lower lower limit used in \code{R} \code{optim} rotuine. Default is \code{c(1e-3,1e-3)}.
#' @param subjects Total number of customers or subject default 1000
#'
#' @return
#' \item{fitted:}{Fitted values based on historical data}
#' \item{projected:}{Projected \code{h} values based on historical data}
#' \item{max.likelihood:}{Maximum Likelihood of Beta Geometric}
#' \item{params - a, b:}{Returns a and b paramters from maximum likelihood estimation for beta distribution}
#'
#' @examples
#' surv_value <- c(100,86.9,74.3,65.3,59.3)
#' h <- 6
#' BG(surv_value,h)
#'
#' @references {Fader P, Hardie B. How to project customer retention. Journal of Interactive Marketing. 2007;21(1):76-90.}
#' @export


BG <- function(surv_value, h, lower = c(1e-3, 1e-3),subjects=1000) {
  surv <- surv_value

  if (surv[1] != 100)
    stop("Starting Value should be 100")

  if (any(surv[-1] >= 100) |
      any(surv[-1] < 0))
    stop("Starting Value should be 100 and non-starting value should be between 0 and less than 100")

  surv <- (surv/100)*subjects

  t <- length(surv)

  die <- diff(-surv)

  s <- rep(NA, length(surv))
  p <- rep(NA, length(surv))

  bg.log.lik <- function(params) {
    a <- params[1]
    b <- params[2]


    i = 0:(t - 1)

    s <- beta(a, b + i) / beta(a, b)

    p <- diff(-s)

    ll_ <- (die[i]) * log(p[i])

    ll <- sum(ll_) + (surv[t]) * log(s[t])

    return(-ll)
  }


  max.lik.sgb <- tryCatch({
    stats::optim(c(1, 2),
                 fn = bg.log.lik,
                 lower = lower,
                 method = "L-BFGS-B")
  }, error = function(error_condition) {
    message(
      "Note: stats::optim not working switching to nloptr::slsqp for maximum likelihood optimization"
    )
    nloptr::slsqp(c(1, 2), fn = bg.log.lik, lower = lower)
  })


  a <- max.lik.sgb$par[1]
  b <- max.lik.sgb$par[2]

  k <- 0:(t+h)

  sbg <- (beta(a, b + k) / beta(a, b)) * 100

  # Vectorized extraction of projected values
  projected <- if (h > 0) {
    projected <- sbg[(t + 1):(t + h)]
  } else {
    message("Forecast horizon h is: ",h,", No Forecast generated.")
    projected <- numeric(0) # Return an empty numeric vector if h = 0
  }

  list(
    fitted = sbg[1:t],
    projected = projected,
    max.likelihood = max.lik.sgb$value,
    params = c(a = a, b = b)
  )


}