equi.pii.LL: P(II) log-likelihood function
In ghuiber/BTYD2: Implementing BTYD Models With The Log Sum Exp Patch

Description Usage Arguments Details Value References Examples

Log-likelihood function for the Pareto distribution of the second kind, P(II), also known as the Lomax distribution.

1	equi.pii.LL(params, x, t, n)

`params`	A vector of the two Gamma parameters r and alpha.
`x`	A vector with the number of transactions at discrete times t.
`t`	The timeline of transactions in vector x.
`n`	The total possible number of transactions at time t.

This is a distribution of continuous-time duration data. Observed duration times follow an exponential distribution with a rate lambda that is specific to each individual; the individual values of the lambda parameter are drawn from a gamma distribution with shape parameter r and inverse-scale parameter alpha. The P(II) pdf and cdf have closed forms where t is conditioned on r and alpha alone.

In other places alpha is called a scale parameter; R documentation calls it a rate parameter.

The log likelihood of observing the combination of transactions x over timeline t given that x follows a P(II) distribution with parameters params.

Peter S. Fader, Bruce G. S. Hardie, Daniel McCarthy & Ramnath Vaidyanathan (2019) Exploring the Equivalence of Two Common Mixture Models for Duration Data, The American Statistician, 73:3, 288-295, DOI: 10.1080/00031305.2018.1543134 Web

# Get the calibration data. Remember that
# counts of households are cumulative, so 
# turn them to incremental:
library(tidyverse)
data("ace_snackfoods")
data <- ace_snackfoods$calibration %>% 
  rename(x = `Number of households`) %>% 
  mutate(bought = x - lag(x)) %>% 
  mutate(bought = if_else(is.na(bought), x, bought))
  
# Get maximum likelihood estimates:   
mlset <- optim(par = c(0, 3), 
               fn = equi.pii.LL, 
               x = data$bought, 
               t = data$Week, 
               n = 1499, 
               method="BFGS", 
               control=list(fnscale=-1))

# Log-likelihood maximum:
LL <- mlset$value

# Parameter estimates:
# Recall that bg.LL exponentiates par in
# order to enforce r > 0, alpha > 0.
params <- exp(mlset$par)