# Copyright (c) 2026 INNOVATIO SAS
# SPDX-License-Identifier: MIT
# MSC-P-046 (marketing-science-center.com): predict how much a customer will spend per purchase.
# Base R only. Same computations, same order and same printed lines as msc-p046-reference.py: the
# gamma-gamma model of spend per transaction (Fader, Hardie and Lee 2005; Fader and Hardie 2013) checked on
# the values the authors publish for the CDNOW sample, then validated on the held-out weeks 40 to 78.
# The data file CDNOW_sample.txt is read from the current directory, next to the program, or from
# editorial/sources/MSC-P-046/ when the program runs inside the repository; it is downloadable from https://www.brucehardie.com/datasets/.

SCRIPT_ARG <- sub("^--file=", "", grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE))
HERE <- if (length(SCRIPT_ARG)) dirname(normalizePath(SCRIPT_ARG[1])) else "."
DATA_CANDIDATES <- c("CDNOW_sample.txt", file.path(HERE, "CDNOW_sample.txt"), file.path(HERE, "..", "..", "sources", "MSC-P-046", "CDNOW_sample.txt"))
RECORDS <- 6919L
CUSTOMERS <- 2357L
CALIBRATION_END <- 272   # days after 1 January 1997: 30 September 1997, end of week 39
HOLDOUT_END <- 545       # 30 June 1998, end of week 78
LCG_SEED <- 20261007
LCG_MODULUS <- 2147483648
LCG_INCREMENT <- 12345
LCG_HIGH <- 16838   # 1103515245 = 16838 * 65536 + 20077, split to keep products exact
LCG_LOW <- 20077
BOOTSTRAP <- 1000L
NM_MAX_ITER <- 20000L
NM_TOL <- 1e-9    # the log-likelihood, a sum of some 3,000 logarithms near -4,000, carries rounding noise of about 1e-11
NM_SIZE <- 1e-7
SIMPSON_STEPS <- 20000L
GRADIENT_STEP <- 1e-4
GRADIENT_TOL <- 1e-3
CV_MIN_PURCHASES <- 3L
CV_SIMULATIONS <- 200L
EXAMPLE_SPEND <- 100.0   # the authors' customer A: one repeat purchase totalling $100
BANDS <- list(c(0.0, 20.0), c(20.0, 50.0), c(50.0, 1e9))

PUB_REPEATERS <- 946L
PUB_TABLE1 <- c("minimum" = 2.99, "25th percentile" = 15.75, "median" = 27.50, "75th percentile" = 41.80,
                "maximum" = 299.63, "mean" = 35.08, "standard deviation" = 30.28, "mode" = 14.96)
PUB_PQG <- c(6.25, 3.74, 15.44)
PUB_MEAN_GAP_CENTS <- 9
PUB_MODEL_MODE <- 19
PUB_OBSERVED_MODE <- 15
PUB_CORRELATION <- 0.11
PUB_CORRELATION_WITHOUT <- 0.06
PUB_P_WITHOUT <- 0.08
PUB_OUTLIER <- c(21, 300)
PUB_LL_39 <- -4659
PUB_LL_39_AT_78 <- -4661
PUB_MEAN_78 <- 36
PUB_SKEWNESS <- 4
PUB_KURTOSIS <- 17

new_stream <- function(seed) {
  env <- new.env()
  env$state <- seed
  env
}
uniform <- function(s) {
  high <- (LCG_HIGH * s$state) %% LCG_MODULUS
  s$state <- (high * 65536 + LCG_LOW * s$state + LCG_INCREMENT) %% LCG_MODULUS
  (s$state + 0.5) / LCG_MODULUS
}

normal <- function(s) {
  u1 <- uniform(s)
  u2 <- uniform(s)
  sqrt(-2.0 * log(u1)) * cos(2.0 * pi * u2)
}

# Gamma(shape, 1) for shape >= 1, Marsaglia and Tsang (2000).
gamma_draw <- function(s, shape) {
  d <- shape - 1.0 / 3.0
  c <- 1.0 / sqrt(9.0 * d)
  repeat {
    x <- normal(s)
    v <- 1.0 + c * x
    if (v <= 0.0) next
    v <- v * v * v
    u <- uniform(s)
    if (log(u) < 0.5 * x * x + d - d * v + d * log(v)) return(d * v)
  }
}

plain_sum <- function(values) {
  total <- 0
  for (v in values) total <- total + v
  total
}

fmt <- function(x, digits) {
  text <- sprintf(paste0("%.", digits, "f"), x)
  if (startsWith(text, "-") && as.numeric(text) == 0) substring(text, 2) else text
}

pct <- function(x, digits = 1) fmt(100.0 * x, digits)

nearest <- function(x) floor(x + 0.5)

fail <- function(message) {
  cat(paste0("error: ", message, "\n"), file = stderr())
  quit(status = 1)
}

out <- function(...) cat(paste0(..., "\n"), sep = "")

# log Gamma written out as in Python; log and exp may differ by one ulp between the two, so the parity holds at the printed precision.
lgam <- function(x) {
  prod <- 1.0
  while (x < 10.0) {
    prod <- prod * x
    x <- x + 1.0
  }
  z <- 1.0 / x
  z2 <- z * z
  series <- z * (1.0 / 12.0 - z2 * (1.0 / 360.0 - z2 * (1.0 / 1260.0 - z2 * (1.0 / 1680.0 - z2 / 1188.0))))
  (x - 0.5) * log(x) - x + 0.5 * log(2.0 * pi) + series - log(prod)
}

betacf <- function(a, b, x) {
  tiny <- 1e-300
  qab <- a + b
  qap <- a + 1.0
  qam <- a - 1.0
  c <- 1.0
  d <- 1.0 - qab * x / qap
  if (abs(d) < tiny) d <- tiny
  d <- 1.0 / d
  h <- d
  for (m in 1:300) {
    m2 <- 2.0 * m
    aa <- m * (b - m) * x / ((qam + m2) * (a + m2))
    d <- 1.0 + aa * d
    if (abs(d) < tiny) d <- tiny
    c <- 1.0 + aa / c
    if (abs(c) < tiny) c <- tiny
    d <- 1.0 / d
    h <- h * (d * c)
    aa <- -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
    d <- 1.0 + aa * d
    if (abs(d) < tiny) d <- tiny
    c <- 1.0 + aa / c
    if (abs(c) < tiny) c <- tiny
    d <- 1.0 / d
    delta <- d * c
    h <- h * delta
    if (abs(delta - 1.0) < 1e-15) break
  }
  h
}

incomplete_beta <- function(a, b, x) {
  if (x <= 0.0) return(0.0)
  if (x >= 1.0) return(1.0)
  front <- exp(lgam(a + b) - lgam(a) - lgam(b) + a * log(x) + b * log(1.0 - x))
  if (x < (a + 1.0) / (a + b + 2.0)) return(front * betacf(a, b, x) / a)
  1.0 - front * betacf(b, a, 1.0 - x) / b
}

t_two_sided <- function(t, df) incomplete_beta(df / 2.0, 0.5, df / (df + t * t))

nelder_mead <- function(f, start, step) {
  k <- length(start)
  pts <- list(start)
  for (j in seq_len(k)) {
    p <- start
    p[j] <- p[j] + step
    pts[[j + 1]] <- p
  }
  vals <- vapply(pts, f, numeric(1))
  converged <- FALSE
  for (iter in seq_len(NM_MAX_ITER)) {
    ord <- order(vals, seq_along(vals))
    pts <- pts[ord]
    vals <- vals[ord]
    size <- 0
    for (i in 2:(k + 1)) for (j in seq_len(k)) size <- max(size, abs(pts[[i]][j] - pts[[1]][j]))
    if (vals[k + 1] - vals[1] < NM_TOL && size < NM_SIZE) {
      converged <- TRUE
      break
    }
    centroid <- numeric(k)
    for (j in seq_len(k)) {
      cc <- 0.0
      for (i in seq_len(k)) cc <- cc + pts[[i]][j]
      centroid[j] <- cc / k
    }
    refl <- centroid + (centroid - pts[[k + 1]])
    fr <- f(refl)
    if (fr < vals[1]) {
      exp_pt <- centroid + 2.0 * (centroid - pts[[k + 1]])
      fe <- f(exp_pt)
      if (fe < fr) {
        pts[[k + 1]] <- exp_pt
        vals[k + 1] <- fe
      } else {
        pts[[k + 1]] <- refl
        vals[k + 1] <- fr
      }
    } else if (fr < vals[k]) {
      pts[[k + 1]] <- refl
      vals[k + 1] <- fr
    } else {
      con <- if (fr < vals[k + 1]) centroid + 0.5 * (refl - centroid) else centroid + 0.5 * (pts[[k + 1]] - centroid)
      fc <- f(con)
      if (fc < min(fr, vals[k + 1])) {
        pts[[k + 1]] <- con
        vals[k + 1] <- fc
      } else {
        for (i in 2:(k + 1)) {
          pts[[i]] <- pts[[1]] + 0.5 * (pts[[i]] - pts[[1]])
          vals[i] <- f(pts[[i]])
        }
      }
    }
  }
  best <- order(vals, seq_along(vals))[1]
  list(point = pts[[best]], value = vals[best], converged = converged)
}

# ---- data -------------------------------------------------------------------------------------------------

day_index <- function(stamp) {
  year <- as.integer(substr(stamp, 1, 4))
  month <- as.integer(substr(stamp, 5, 6))
  day <- as.integer(substr(stamp, 7, 8))
  before <- c(0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334)
  365 * (year - 1997) + before[month] + day - 1
}

read_data <- function() {
  found <- DATA_CANDIDATES[file.exists(DATA_CANDIDATES)]
  if (!length(found)) fail("CDNOW_sample.txt not found: download it from https://www.brucehardie.com/datasets/")
  lines <- readLines(found[1], warn = FALSE)
  lines <- lines[nzchar(trimws(lines))]
  fields <- strsplit(trimws(lines), "[[:space:]]+")
  data.frame(
    cid = vapply(fields, function(f) as.integer(f[2]), integer(1)),
    day = vapply(fields, function(f) day_index(f[3]), numeric(1)),
    dollars = vapply(fields, function(f) as.numeric(f[5]), numeric(1))
  )
}

purchase_days <- function(records) {
  days <- vector("list", CUSTOMERS)
  for (i in seq_len(nrow(records))) {
    cid <- records$cid[i]
    cust <- days[[cid]]
    n <- if (is.null(cust)) 0L else nrow(cust)
    if (n > 0L && cust[n, 1] == records$day[i]) {
      cust[n, 2] <- cust[n, 2] + records$dollars[i]
    } else {
      cust <- rbind(cust, c(records$day[i], records$dollars[i]))
    }
    days[[cid]] <- cust
  }
  days
}

summary_all <- function(days, end) {
  x <- integer(CUSTOMERS)
  z <- numeric(CUSTOMERS)
  for (cid in seq_len(CUSTOMERS)) {
    cust <- days[[cid]]
    values <- cust[cust[, 1] <= end, 2]
    x[cid] <- length(values)
    z[cid] <- plain_sum(values) / length(values)
  }
  list(x = x, z = z)
}

summary_rows <- function(days, start, end) {
  x <- integer(CUSTOMERS)
  z <- numeric(CUSTOMERS)
  for (cid in seq_len(CUSTOMERS)) {
    cust <- days[[cid]]
    if (nrow(cust) > 1L) {
      later <- cust[-1, , drop = FALSE]
      values <- later[later[, 1] >= start & later[, 1] <= end, 2]
    } else {
      values <- numeric(0)
    }
    x[cid] <- length(values)
    z[cid] <- if (length(values)) plain_sum(values) / length(values) else 0.0
  }
  list(x = x, z = z)
}

# ---- the gamma-gamma model --------------------------------------------------------------------------------

repeaters <- function(x, z) {
  distinct <- sort(unique(x))
  groups <- lapply(distinct, function(k) {
    zs <- z[x == k]
    list(k = k, count = as.numeric(length(zs)), sum_lz = plain_sum(log(zs)), zs = zs, lx = log(k))
  })
  list(x = x, z = z, distinct = distinct, groups = groups)
}

loglik <- function(theta, data) {
  p <- theta[1]
  q <- theta[2]
  g <- theta[3]
  lq <- lgam(q)
  lng <- log(g)
  total <- 0.0
  for (grp in data$groups) {
    k <- grp$k
    logs <- plain_sum(log(g + k * grp$zs))
    a <- lgam(p * k + q) - lgam(p * k) - lq + q * lng + p * k * grp$lx
    total <- total + (grp$count * a + (p * k - 1.0) * grp$sum_lz - (p * k + q) * logs)
  }
  total
}

fit <- function(data, start, step) {
  objective <- function(v) -loglik(exp(v), data)
  run <- nelder_mead(objective, start, step)
  point <- run$point
  value <- run$value
  for (i in 1:10) {
    nxt <- nelder_mead(objective, point, 0.1)
    stable <- value - nxt$value < 1e-9
    point <- nxt$point
    value <- nxt$value
    if (stable && nxt$converged) return(list(theta = exp(point), ll = -value))
  }
  fail("the simplex does not stabilise")
}

true_mean_sd <- function(theta) theta[1] * theta[3] / ((theta[2] - 1.0) * sqrt(theta[2] - 2.0))

population_mean <- function(theta) theta[1] * theta[3] / (theta[2] - 1.0)

own_weight <- function(theta, x) theta[1] * x / (theta[1] * x + theta[2] - 1.0)

conditional_mean <- function(theta, x, zbar) {
  if (x == 0) return(population_mean(theta))
  theta[1] * (theta[3] + x * zbar) / (theta[1] * x + theta[2] - 1.0)
}

density_zbar <- function(theta, x, zbar) {
  p <- theta[1]
  q <- theta[2]
  g <- theta[3]
  exp(-lgam(p * x) - lgam(q) + lgam(p * x + q) + q * log(g / (g + x * zbar)) + p * x * log(x * zbar / (g + x * zbar))) / zbar
}

model_mode <- function(theta, counts) {
  ks <- as.integer(names(counts))
  total <- plain_sum(as.numeric(counts))
  best <- -1.0
  best_y <- 0
  for (y in 1:300) {
    terms <- numeric(length(ks))
    for (i in seq_along(ks)) terms[i] <- counts[[i]] * density_zbar(theta, ks[i], as.numeric(y))
    f <- plain_sum(terms) / total
    if (f > best) {
      best <- f
      best_y <- y
    }
  }
  best_y
}

# ---- descriptive statistics -------------------------------------------------------------------------------

mean_abs <- function(rows, key) plain_sum(abs(rows$z2 - rows[[key]])) / nrow(rows)

median_of <- function(xs) {
  s <- sort(xs)
  k <- length(s)
  if (k %% 2 == 1) s[k %/% 2 + 1] else (s[k %/% 2] + s[k %/% 2 + 1]) / 2.0
}

mean_sd <- function(xs) {
  m <- plain_sum(xs) / length(xs)
  c(m, sqrt(plain_sum((xs - m) * (xs - m)) / (length(xs) - 1)))
}

weibull_quantile <- function(sorted_xs, prob) {
  h <- (length(sorted_xs) + 1) * prob
  lo <- floor(h)
  sorted_xs[lo] + (h - lo) * (sorted_xs[lo + 1] - sorted_xs[lo])
}

moments <- function(xs) {
  m <- plain_sum(xs) / length(xs)
  d <- xs - m
  m2 <- plain_sum(d * d) / length(xs)
  m3 <- plain_sum(d * d * d) / length(xs)
  m4 <- plain_sum((d * d) * (d * d)) / length(xs)
  c(m3 / (m2 * sqrt(m2)), m4 / (m2 * m2) - 3.0)
}

cents_mode <- function(xs) {
  keys <- floor(xs * 100.0 + 0.5)
  tab <- table(keys)
  top <- max(tab)
  c(min(as.numeric(names(tab)[tab == top])) / 100.0, top)
}

correlation <- function(xs, ys) {
  mx <- plain_sum(xs) / length(xs)
  my <- plain_sum(ys) / length(ys)
  sxy <- plain_sum((xs - mx) * (ys - my))
  sxx <- plain_sum((xs - mx) * (xs - mx))
  syy <- plain_sum((ys - my) * (ys - my))
  sxy / sqrt(sxx * syy)
}

quantile_pair <- function(xs) {
  s <- sort(xs)
  c(s[nearest(0.025 * length(s))], s[nearest(0.975 * length(s))])
}

simpson <- function(f, a, b, steps) {
  h <- (b - a) / steps
  total <- f(a) + f(b)
  for (i in 1:(steps - 1)) total <- total + (if (i %% 2 == 1) 4.0 else 2.0) * f(a + i * h)
  total * h / 3.0
}

# ---- checks ---------------------------------------------------------------------------------------------------

check_invariants <- function(theta) {
  checks <- 0
  ensure <- function(condition, message) {
    checks <<- checks + 1
    if (!condition) fail(paste0("invariant failed: ", message))
  }
  ensure(abs(lgam(1.0)) < 1e-13 && abs(lgam(2.0)) < 1e-13, "log Gamma(1) = log Gamma(2) = 0")
  ensure(abs(lgam(0.5) - 0.5 * log(pi)) < 1e-13, "log Gamma(1/2) = log(pi) / 2")
  for (v in c(0.3, 3.7441, 6.2498, 43.75, 131.25)) {
    ensure(abs(lgam(v) - lgamma(v)) < 1e-11 * max(1.0, abs(lgamma(v))), paste0("declared log Gamma(", v, ") against the library"))
  }
  ensure(abs(lgam(7.3) - lgam(6.3) - log(6.3)) < 1e-12, "log Gamma(x + 1) = log Gamma(x) + log x")
  ensure(abs(t_two_sided(1.959963984540054, 1e7) - 0.05) < 1e-6, "two-sided p of 1.96 with many degrees of freedom is 0.05")
  ensure(abs(incomplete_beta(2.0, 3.0, 0.4) - 0.5248) < 1e-12, "I_0.4(2, 3) = 0.5248")
  p <- theta[1]
  q <- theta[2]
  g <- theta[3]
  for (x in c(1, 4)) {
    a <- p * x
    b <- q
    norm <- lgam(a + b) - lgam(a) - lgam(b)
    dens <- function(u) if (u <= 0.0 || u >= 1.0) 0.0 else exp(norm + (a - 1.0) * log(u) + (b - 1.0) * log(1.0 - u))
    ensure(abs(simpson(dens, 0.0, 1.0, SIMPSON_STEPS) - 1.0) < 1e-8, paste0("f(zbar | x = ", x, ") integrates to 1"))
    mean <- simpson(function(u) if (u < 1.0) dens(u) * u / (1.0 - u) * g / x else 0.0, 0.0, 1.0, SIMPSON_STEPS)
    ensure(abs(mean - population_mean(theta)) < 1e-4, paste0("E(zbar | x = ", x, ") equals p gamma / (q - 1)"))
    zbar <- 27.5
    form_a <- exp(lgam(a + q) - lgam(a) - lgam(q) + (a - 1.0) * log(zbar) + a * log(x) + q * log(g) - (a + q) * log(g + x * zbar))
    ensure(abs(form_a - density_zbar(theta, x, zbar)) < 1e-12 * form_a, paste0("forms (1a) and (1b) agree at x = ", x))
    w <- own_weight(theta, x)
    ensure(abs(conditional_mean(theta, x, zbar) - ((1.0 - w) * population_mean(theta) + w * zbar)) < 1e-10, paste0("equation (5) is a weighted average at x = ", x))
  }
  ensure(own_weight(theta, 1000000) > 0.9999, "the weight of the own average tends to 1")
  checks
}

# Stops unless the fitted point is a stationary point: the log-likelihood gradient on log p, log q, log gamma is near zero.
check_gradient <- function(theta, data) {
  point <- log(theta)
  for (j in 1:3) {
    up <- point
    down <- point
    up[j] <- up[j] + GRADIENT_STEP
    down[j] <- down[j] - GRADIENT_STEP
    slope <- (loglik(exp(up), data) - loglik(exp(down), data)) / (2.0 * GRADIENT_STEP)
    if (abs(slope) > GRADIENT_TOL) fail("the fitted point is not a maximum of the log-likelihood")
  }
}

published_checks <- function(n_customers, rdata, theta, ll, theta78, ll39_at_78, counts) {
  z <- rdata$z
  s <- sort(z)
  msd <- mean_sd(z)
  mode <- cents_mode(z)[1]
  table1 <- c("minimum" = s[1], "25th percentile" = weibull_quantile(s, 0.25), "median" = weibull_quantile(s, 0.5),
              "75th percentile" = weibull_quantile(s, 0.75), "maximum" = s[length(s)], "mean" = msd[1], "standard deviation" = msd[2], "mode" = mode)
  reproduced <- 0
  out("Published values, Fader and Hardie (2013) and Fader, Hardie and Lee (2005), on the same data:")
  out("  repeat buyers in weeks 1-39: ", length(z), " of ", n_customers, " (published ", PUB_REPEATERS, " of ", CUSTOMERS, ")")
  if (length(z) != PUB_REPEATERS || n_customers != CUSTOMERS) fail("the number of repeat buyers is not reproduced")
  for (name in names(table1)) {
    ok <- abs(table1[[name]] - PUB_TABLE1[[name]]) <= 0.0051
    reproduced <- reproduced + ok
    out("  Table 1, ", name, ": ", fmt(table1[[name]], 4), " (published ", fmt(PUB_TABLE1[[name]], 2), ") ", if (ok) "reproduced" else "NOT reproduced")
  }
  if (reproduced != length(table1)) fail("Table 1 is not reproduced")
  mom <- moments(z)
  out("  2005 Table 1, skewness: ", fmt(mom[1], 2), " (published ", PUB_SKEWNESS, ") ", if (nearest(mom[1]) == PUB_SKEWNESS) "reproduced" else "NOT reproduced")
  out("  2005 Table 1, excess kurtosis: ", fmt(mom[2], 2), " (published ", PUB_KURTOSIS, ") ", if (nearest(mom[2]) == PUB_KURTOSIS) "reproduced" else "NOT reproduced")
  if (nearest(mom[2]) != PUB_KURTOSIS) fail("the kurtosis is not reproduced")
  names_pqg <- c("p", "q", "gamma")
  for (i in 1:3) {
    out("  estimate ", names_pqg[i], ": ", fmt(theta[i], 4), " (published ", fmt(PUB_PQG[i], 2), ")")
    if (abs(theta[i] - PUB_PQG[i]) > 0.0051) fail(paste0("the estimate of ", names_pqg[i], " is not reproduced"))
  }
  gap <- population_mean(theta) - msd[1]
  out("  theoretical mean ", fmt(population_mean(theta), 4), " minus observed mean ", fmt(msd[1], 4), ": ", fmt(gap, 4), " dollar (published: ", PUB_MEAN_GAP_CENTS, " cents)")
  if (nearest(100.0 * gap) != PUB_MEAN_GAP_CENTS) fail("the gap between theoretical and observed mean is not reproduced")
  model <- model_mode(theta, counts)
  out("  mode of the model density of zbar: ", model, " dollars (published ", PUB_MODEL_MODE, "); observed mode ", fmt(mode, 2), " (published ", PUB_OBSERVED_MODE, ")")
  if (model != PUB_MODEL_MODE || nearest(mode) != PUB_OBSERVED_MODE) fail("the modes are not reproduced")
  xs <- as.numeric(rdata$x)
  r <- correlation(xs, z)
  outlier <- which(rdata$x == PUB_OUTLIER[1])
  if (length(outlier) != 1 || nearest(z[outlier]) != PUB_OUTLIER[2]) fail("the outlier (21 transactions, $300) is not found")
  keep <- seq_along(z)[-outlier]
  r2 <- correlation(xs[keep], z[keep])
  df <- length(keep) - 2
  pvalue <- t_two_sided(r2 * sqrt(df) / sqrt(1.0 - r2 * r2), as.numeric(df))
  out("  correlation of frequency and mean spend: ", fmt(r, 4), " (published ", fmt(PUB_CORRELATION, 2), "); without the customer with ", rdata$x[outlier], " purchases and a mean of ", fmt(z[outlier], 2), ": ", fmt(r2, 4), ", p = ", fmt(pvalue, 4), " (published ", fmt(PUB_CORRELATION_WITHOUT, 2), ", p = ", fmt(PUB_P_WITHOUT, 2), ")")
  if (nearest(100.0 * r) != nearest(100.0 * PUB_CORRELATION) || nearest(100.0 * r2) != nearest(100.0 * PUB_CORRELATION_WITHOUT) || nearest(100.0 * pvalue) != nearest(100.0 * PUB_P_WITHOUT)) fail("the correlations are not reproduced")
  sum_log_x <- plain_sum(log(rdata$x))
  out("  log-likelihood of the mean spends, weeks 1-39: ", fmt(ll, 4), "; minus the sum of log x (", fmt(sum_log_x, 4), "): ", fmt(ll - sum_log_x, 4), " (published ", PUB_LL_39, ")")
  out("  78-week estimates p = ", fmt(theta78[1], 4), ", q = ", fmt(theta78[2], 4), ", gamma = ", fmt(theta78[3], 4), "; their weeks 1-39 log-likelihood: ", fmt(ll39_at_78, 4), "; minus the sum of log x: ", fmt(ll39_at_78 - sum_log_x, 4), " (published ", PUB_LL_39_AT_78, ")")
  out("  78-week population mean: ", fmt(population_mean(theta78), 4), " (published about ", PUB_MEAN_78, ")")
  if (nearest(ll - sum_log_x) != PUB_LL_39 || nearest(ll39_at_78 - sum_log_x) != PUB_LL_39_AT_78 || nearest(population_mean(theta78)) != PUB_MEAN_78) fail("the 2005 log-likelihoods or the 78-week mean are not reproduced")
  sum_log_x
}

# ---- main -------------------------------------------------------------------------------------------------

main <- function() {
  records <- read_data()
  customers <- sort(unique(records$cid))
  if (nrow(records) != RECORDS || !identical(customers, seq_len(CUSTOMERS))) fail("CDNOW_sample.txt is not the published file (6,919 records, customers 1 to 2,357)")
  days <- purchase_days(records)
  if (max(vapply(days, function(d) d[1, 1], numeric(1))) > 89) fail("a customer's first purchase falls after March 1997")
  purchases <- sum(vapply(days, nrow, integer(1)))
  dollars <- plain_sum(records$dollars)
  out("MSC-P-046: gamma-gamma model of spend per purchase, CDNOW sample (Fader and Hardie)")
  out("Data: ", nrow(records), " records, ", CUSTOMERS, " customers, ", purchases, " purchase days once same-day transactions are added up, ", fmt(dollars, 2), " dollars")
  out("Calibration: weeks 1-39 (to 30 September 1997); holdout: weeks 40-78 (1 October 1997 to 30 June 1998)")

  cal <- summary_rows(days, 0, CALIBRATION_END)
  full <- summary_rows(days, 0, HOLDOUT_END)
  hold <- summary_rows(days, CALIBRATION_END + 1, HOLDOUT_END)
  rdata <- repeaters(cal$x[cal$x > 0], cal$z[cal$x > 0])
  counts <- table(rdata$x)
  fitted <- fit(rdata, c(0.0, 0.0, 0.0), 1.0)
  theta <- fitted$theta
  ll <- fitted$ll
  full_rep <- repeaters(full$x[full$x > 0], full$z[full$x > 0])
  theta78 <- fit(full_rep, c(0.0, 0.0, 0.0), 1.0)$theta
  ll39_at_78 <- loglik(theta78, rdata)

  check_gradient(theta, rdata)
  published_checks(length(cal$x), rdata, theta, ll, theta78, ll39_at_78, counts)
  out("Invariants: ", check_invariants(theta), " checks passed")

  p <- theta[1]
  q <- theta[2]
  g <- theta[3]
  ez <- population_mean(theta)
  out("")
  out("Fitted on weeks 1-39 (946 repeat buyers):")
  out("  p = ", fmt(p, 4), ", q = ", fmt(q, 4), ", gamma = ", fmt(g, 4), ", log-likelihood ", fmt(ll, 2))
  out("  coefficient of variation of a customer's purchases around their own mean, 1 / sqrt(p): ", fmt(1.0 / sqrt(p), 4))
  out("  population mean spend per purchase E(Z) = p gamma / (q - 1): ", fmt(ez, 4), " dollars")
  out("  spread of the customers' true means, standard deviation: ", fmt(true_mean_sd(theta), 4), " dollars")
  out("  repeat buyers by number of repeat purchases: ", paste0(names(counts), ": ", as.integer(counts), collapse = ", "))

  out("")
  out("Weight of the customer's own average in E(Z | x, zbar), and examples:")
  for (x in c(1, 2, 3, 4, 6, 8, 10)) {
    out("  x = ", x, ": weight ", fmt(own_weight(theta, x), 4), " (78-week fit ", fmt(own_weight(theta78, x), 4), "); zbar 20 -> ", fmt(conditional_mean(theta, x, 20.0), 2), ", zbar 50 -> ", fmt(conditional_mean(theta, x, 50.0), 2), ", zbar 100 -> ", fmt(conditional_mean(theta, x, 100.0), 2))
  }
  first90 <- which(vapply(1:99, function(x) own_weight(theta, x) >= 0.9, logical(1)))[1]
  first90_78 <- which(vapply(1:99, function(x) own_weight(theta78, x) >= 0.9, logical(1)))[1]
  out("  repeat purchases before the own average weighs 90 %: ", first90, " (78-week fit: ", first90_78, ")")
  out("  customer A, one repeat purchase of ", fmt(EXAMPLE_SPEND, 0), " dollars: E(Z) = ", fmt(conditional_mean(theta, 1, EXAMPLE_SPEND), 2), " dollars")
  out("  customer A with the 78-week fit: E(Z) = ", fmt(conditional_mean(theta78, 1, EXAMPLE_SPEND), 2), " dollars")

  alls <- summary_all(days, CALIBRATION_END)
  keep_all <- alls$z > 0.0
  all_fit_set <- repeaters(alls$x[keep_all], alls$z[keep_all])
  fitted_all <- fit(all_fit_set, c(0.0, 0.0, 0.0), 1.0)
  theta_all <- fitted_all$theta
  out("")
  out("Variant of the page: gamma-gamma fitted on all purchases of weeks 1-39, the first one included:")
  out("  customers with a mean of zero left out of the fit: ", CUSTOMERS - sum(keep_all))
  out("  p = ", fmt(theta_all[1], 4), ", q = ", fmt(theta_all[2], 4), ", gamma = ", fmt(theta_all[3], 4), ", log-likelihood ", fmt(fitted_all$ll, 2), "; population mean ", fmt(population_mean(theta_all), 4), " dollars")
  first_values <- vapply(days, function(d) d[1, 2], numeric(1))
  repeat_values <- unlist(lapply(days, function(d) if (nrow(d) > 1L) { later <- d[-1, , drop = FALSE]; later[later[, 1] <= CALIBRATION_END, 2] } else numeric(0)))
  holdout_values <- unlist(lapply(days, function(d) d[d[, 1] > CALIBRATION_END & d[, 1] <= HOLDOUT_END, 2]))
  out("  mean value of a first purchase ", fmt(plain_sum(first_values) / length(first_values), 2), ", of a repeat purchase in weeks 1-39 ", fmt(plain_sum(repeat_values) / length(repeat_values), 2), ", of a purchase in weeks 40-78 ", fmt(plain_sum(holdout_values) / length(holdout_values), 2), " (", length(holdout_values), " purchases)")
  first_rep <- first_values[cal$x > 0]
  first_none <- first_values[cal$x == 0]
  ratios <- first_values[cal$x > 0] / cal$z[cal$x > 0]
  out("  first purchase of the ", length(first_rep), " repeat buyers ", fmt(plain_sum(first_rep) / length(first_rep), 2), " against the mean of their repeat purchases ", fmt(plain_sum(rdata$z) / length(rdata$z), 2), " (median ratio ", fmt(median_of(ratios), 4), "); first purchase of the ", length(first_none), " customers without repeat purchase ", fmt(plain_sum(first_none) / length(first_none), 2))
  ms_rep <- mean_sd(first_rep)
  ms_none <- mean_sd(first_none)
  se_gap <- sqrt(ms_rep[2] * ms_rep[2] / length(first_rep) + ms_none[2] * ms_none[2] / length(first_none))
  out("  first purchase, repeat buyers minus customers without repeat purchase: ", fmt(ms_rep[1] - ms_none[1], 2), " (standard error ", fmt(se_gap, 2), ")")
  out("  one purchase of ", fmt(EXAMPLE_SPEND, 0), " dollars: ", fmt(conditional_mean(theta_all, 1, EXAMPLE_SPEND), 2), "; two purchases averaging ", fmt(EXAMPLE_SPEND, 0), ": ", fmt(conditional_mean(theta_all, 2, EXAMPLE_SPEND), 2), "; weight of one purchase ", fmt(own_weight(theta_all, 1), 4), ", of two ", fmt(own_weight(theta_all, 2), 4))

  out("")
  out("Holdout, weeks 40-78 (error = actual mean spend per purchase minus prediction; standard errors across customers):")
  rows <- data.frame(x = cal$x, zbar = cal$z, n_all = alls$x, z_all = alls$z, n2 = hold$x, z2 = hold$z)
  rows$rep <- rows$zbar
  rows$all <- rows$z_all
  rows$pop <- ez
  rows$gg <- vapply(seq_len(CUSTOMERS), function(i) conditional_mean(theta, rows$x[i], rows$zbar[i]), numeric(1))
  rows$gg_all <- vapply(seq_len(CUSTOMERS), function(i) conditional_mean(theta_all, rows$n_all[i], rows$z_all[i]), numeric(1))
  labels <- c(rep = "average of repeat purchases", all = "average of all purchases", pop = "population mean",
              gg = "gamma-gamma, repeat purchases (authors)", gg_all = "gamma-gamma, all purchases (page)")
  sets <- list(
    list(label = "repeat buyers of weeks 1-39 who buy again", test = rows[rows$x > 0 & rows$n2 > 0, ], keys = c("rep", "all", "pop", "gg", "gg_all")),
    list(label = "customers without repeat purchase in weeks 1-39 who buy", test = rows[rows$x == 0 & rows$n2 > 0, ], keys = c("all", "pop", "gg_all")),
    list(label = "all customers who buy in weeks 40-78", test = rows[rows$n2 > 0, ], keys = c("all", "gg", "gg_all"))
  )
  for (st in sets) {
    test <- st$test
    out("  ", st$label, ": ", nrow(test))
    for (key in st$keys) {
      e <- test$z2 - test[[key]]
      mae <- plain_sum(abs(e)) / length(e)
      rmse <- sqrt(plain_sum(e * e) / length(e))
      ms <- mean_sd(e)
      weighted <- plain_sum(test$n2 * abs(test$z2 - test[[key]])) / plain_sum(as.numeric(test$n2))
      line <- paste0("    ", labels[[key]], ": mean absolute error ", fmt(mae, 2), ", root mean square error ", fmt(rmse, 2), ", mean error ", fmt(ms[1], 2), " (standard error ", fmt(ms[2] / sqrt(length(e)), 2), "), mean absolute error weighted by purchases ", fmt(weighted, 2))
      if (key != "all") {
        d <- abs(test$z2 - test[[key]]) - abs(test$z2 - test$all)
        dms <- mean_sd(d)
        line <- paste0(line, "; minus average of all purchases ", fmt(dms[1], 2), " (standard error ", fmt(dms[2] / sqrt(length(d)), 2), ")")
      }
      if (!(key %in% c("gg", "pop"))) {
        d <- abs(test$z2 - test[[key]]) - abs(test$z2 - test$gg)
        dms <- mean_sd(d)
        line <- paste0(line, "; minus authors' model ", fmt(dms[1], 2), " (standard error ", fmt(dms[2] / sqrt(length(d)), 2), ")")
      }
      out(line)
    }
  }
  test <- sets[[1]]$test
  better <- sum(abs(test$z2 - test$gg_all) < abs(test$z2 - test$all))
  better_rep <- sum(abs(test$z2 - test$gg) < abs(test$z2 - test$rep))
  out("  repeat buyers for whom the page's variant beats the average of all purchases: ", better, " of ", nrow(test), " (", pct(better / nrow(test)), " %); the authors' model beats the average of repeat purchases for ", better_rep, " (", pct(better_rep / nrow(test)), " %)")
  out("  authors' model against average of repeat purchases: mean absolute error ", pct(mean_abs(test, "gg") / mean_abs(test, "rep") - 1.0), " %; page's variant against average of all purchases: ", pct(mean_abs(test, "gg_all") / mean_abs(test, "all") - 1.0), " %")
  out("  by average of repeat purchases in weeks 1-39 (customers, mean zbar, mean prediction of the authors' model, of the page's variant, actual mean spend in weeks 40-78, actual minus authors' prediction with its standard error):")
  for (band in BANDS) {
    sel <- test[test$zbar >= band[1] & test$zbar < band[2], ]
    label <- if (band[2] < 1e8) paste0(fmt(band[1], 0), " to ", fmt(band[2], 0)) else paste0(fmt(band[1], 0), " and more")
    dms <- mean_sd(sel$z2 - sel$gg)
    out("    ", label, ": ", nrow(sel), ", ", fmt(plain_sum(sel$zbar) / nrow(sel), 2), ", ", fmt(plain_sum(sel$gg) / nrow(sel), 2), ", ", fmt(plain_sum(sel$gg_all) / nrow(sel), 2), ", ", fmt(plain_sum(sel$z2) / nrow(sel), 2), ", ", fmt(dms[1], 2), " (", fmt(dms[2] / sqrt(nrow(sel)), 2), ")")
  }
  out("  by repeat purchases in weeks 1-39, mean total spend in weeks 40-78 per customer (customers, actual, authors' model times actual purchases, difference with its standard error):")
  for (lo in 0:7) {
    hi <- if (lo < 7) lo else 1000000
    grp <- rows[rows$x >= lo & rows$x <= hi, ]
    actual <- plain_sum(grp$n2 * grp$z2) / nrow(grp)
    expected <- plain_sum(grp$n2 * grp$gg) / nrow(grp)
    dms <- mean_sd(grp$n2 * (grp$z2 - grp$gg))
    out("    ", if (lo == hi) lo else paste0(lo, "+"), ": ", nrow(grp), ", ", fmt(actual, 2), ", ", fmt(expected, 2), ", ", fmt(dms[1], 2), " (", fmt(dms[2] / sqrt(nrow(grp)), 2), ")")
  }
  actual_total <- plain_sum(rows$n2 * rows$z2)
  for (key in c("gg", "gg_all", "all")) {
    total <- plain_sum(rows$n2 * rows[[key]])
    out("  total holdout spend, ", labels[[key]], " times actual purchases: ", fmt(total, 2), " against ", fmt(actual_total, 2), " (", pct(total / actual_total - 1.0), " %)")
  }

  stream_cv <- new_stream(LCG_SEED + 1)
  cv_k <- integer(0)
  cv_m <- numeric(0)
  cv_v <- numeric(0)
  for (cid in seq_len(CUSTOMERS)) {
    d <- days[[cid]]
    values <- if (nrow(d) > 1L) { later <- d[-1, , drop = FALSE]; later[later[, 1] <= CALIBRATION_END, 2] } else numeric(0)
    if (length(values) >= CV_MIN_PURCHASES) {
      ms <- mean_sd(values)
      cv_k <- c(cv_k, length(values))
      cv_m <- c(cv_m, ms[1])
      cv_v <- c(cv_v, ms[2] / ms[1])
    }
  }
  observed <- median_of(cv_v)
  medians <- numeric(CV_SIMULATIONS)
  for (s in seq_len(CV_SIMULATIONS)) {
    sim <- numeric(length(cv_k))
    for (i in seq_along(cv_k)) {
      draw <- numeric(cv_k[i])
      for (j in seq_len(cv_k[i])) draw[j] <- gamma_draw(stream_cv, p)
      ms <- mean_sd(draw)
      sim[i] <- ms[2] / ms[1]
    }
    medians[s] <- median_of(sim)
  }
  qp <- quantile_pair(medians)
  out("")
  out("Spread of a customer's purchases around their own mean, repeat buyers with at least ", CV_MIN_PURCHASES, " repeat purchases in weeks 1-39 (", length(cv_k), " customers):")
  out("  median observed coefficient of variation ", fmt(observed, 4), "; expected under the model with the same numbers of purchases, mean of ", CV_SIMULATIONS, " simulations ", fmt(plain_sum(medians) / length(medians), 4), ", 95 % of simulations between ", fmt(qp[1], 4), " and ", fmt(qp[2], 4))
  for (band in BANDS) {
    sel <- cv_v[cv_m >= band[1] & cv_m < band[2]]
    label <- if (band[2] < 1e8) paste0(fmt(band[1], 0), " to ", fmt(band[2], 0)) else paste0(fmt(band[1], 0), " and more")
    out("    mean spend ", label, ": ", length(sel), " customers, median coefficient of variation ", fmt(median_of(sel), 4))
  }

  stream <- new_stream(LCG_SEED)
  n <- length(rdata$x)
  names_draws <- c("p", "q", "gamma", "E(Z)", "customer A", "weight x = 1", "coefficient of variation 1 / sqrt(p)", "standard deviation of the true means")
  draws <- matrix(0, nrow = BOOTSTRAP, ncol = length(names_draws))
  start <- log(theta)
  for (b in seq_len(BOOTSTRAP)) {
    idx <- integer(n)
    for (i in seq_len(n)) idx[i] <- floor(uniform(stream) * n) + 1
    th <- fit(repeaters(rdata$x[idx], rdata$z[idx]), start, 0.1)$theta
    draws[b, ] <- c(th[1], th[2], th[3], population_mean(th), conditional_mean(th, 1, EXAMPLE_SPEND), own_weight(th, 1), 1.0 / sqrt(th[1]), if (th[2] > 2.0) true_mean_sd(th) else 1e9)
  }
  out("")
  out("Bootstrap, ", BOOTSTRAP, " resamples of the ", n, " repeat buyers, authors' model refitted (seed ", sprintf("%d", LCG_SEED), "), nominal 95 % interval:")
  for (j in seq_along(names_draws)) {
    qp <- quantile_pair(draws[, j])
    out("  ", names_draws[j], ": ", fmt(qp[1], 4), " to ", fmt(qp[2], 4))
  }
}

main()