# Copyright (c) 2026 INNOVATIO SAS
# SPDX-License-Identifier: MIT
# MSC-P-044 (marketing-science-center.com): value customers who can only buy at fixed occasions, with BG/BB.
# Base R only. Same declared recursion, design and algorithm as msc-p044-reference.py; the two programs
# must print exactly the same lines. Floating-point sums are written as explicit left-to-right loops,
# never with sum(), which accumulates in long double; sum() is used for integer counts only.
LCG_SEED <- 20261005
LCG_MODULUS <- 2147483648
LCG_INCREMENT <- 12345
LCG_HIGH <- 16838   # 1103515245 = 16838 * 65536 + 20077, split to keep products exact
LCG_LOW <- 20077

CUSTOMERS <- 4000L
TRUE_ALPHA <- 1.0
TRUE_BETA <- 0.8
TRUE_GAMMA <- 0.5
TRUE_DELTA <- 2.5
CALIBRATION <- 6L
HOLDOUT <- 5L
MARGIN <- 150
DISCOUNT <- 0.10
TERMS <- 2000L
BOOTSTRAP <- 1000L
REPLICATIONS <- 100L
NM_MAX_ITER <- 20000L
NM_TOL <- 1e-10
# Far from any optimum (a parameter above e^50), the log-likelihood is not evaluated: the simplex is told
# "much worse" instead of meeting an overflow. The same bound is used in Python.
LOG_PARAMETER_MAX <- 50
PENALTY <- 1e300

# Fader, Hardie and Shang (2010), Table 2: 1995 cohort, n = 6; (x, t_x, number of donors).
DONORS <- data.frame(
  x = c(6, 5, 4, 3, 2, 1, 5, 4, 3, 2, 1, 4, 3, 2, 1, 3, 2, 1, 2, 1, 1, 0),
  tx = c(6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 4, 4, 4, 4, 3, 3, 3, 2, 2, 1, 0),
  f = c(1203, 728, 512, 357, 234, 129, 335, 284, 225, 173, 119, 240, 181, 155, 78, 322, 255, 129, 613, 277, 1091, 3464)
)

lcg_state <- LCG_SEED
uniform <- function() {
  high <- (LCG_HIGH * lcg_state) %% LCG_MODULUS
  lcg_state <<- (high * 65536 + LCG_LOW * lcg_state + LCG_INCREMENT) %% LCG_MODULUS
  (lcg_state + 0.5) / LCG_MODULUS
}

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

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

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

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

# The n(n + 1)/2 + 1 recency/frequency patterns (x, t_x), in the order of Table 2.
patterns <- function(n) {
  x <- integer(0)
  tx <- integer(0)
  for (t in n:1) for (k in t:1) {
    x <- c(x, k)
    tx <- c(tx, t)
  }
  data.frame(x = c(x, 0L), tx = c(tx, 0L))
}

# B(a + x, b + m - x) / B(a, b): x purchases in m opportunities, by ascending products.
ratio_ab <- function(a, b, x, m) {
  num <- 1
  if (x > 0) for (k in 0:(x - 1)) num <- num * (a + k)
  if (m - x > 0) for (k in 0:(m - x - 1)) num <- num * (b + k)
  den <- 1
  if (m > 0) for (k in 0:(m - 1)) den <- den * (a + b + k)
  num / den
}

# B(g, d + t) / B(g, d): still alive after t opportunities.
alive_through <- function(g, d, t) {
  s <- 1
  if (t > 0) for (k in 0:(t - 1)) s <- s * ((d + k) / (g + d + k))
  s
}

# B(g + 1, d + t) / B(g, d): alive through t opportunities, dead at the start of opportunity t + 1.
dies_at <- function(g, d, t) alive_through(g, d, t) * g / (g + d + t)

# Equation (5) of the paper, written with ascending products instead of beta functions.
likelihood <- function(par, x, tx, n) {
  total <- ratio_ab(par[1], par[2], x, n) * alive_through(par[3], par[4], n)
  if (n - tx > 0) for (i in 0:(n - tx - 1)) {
    total <- total + ratio_ab(par[1], par[2], x, tx + i) * dies_at(par[3], par[4], tx + i)
  }
  total
}

loglik <- function(par, data, n) {
  total <- 0
  for (j in seq_len(nrow(data))) {
    if (data$f[j] > 0) total <- total + data$f[j] * log(likelihood(par, data$x[j], data$tx[j], n))
  }
  total
}

# Beta-Bernoulli: no dropout, every customer alive for ever.
bb_loglik <- function(a, b, data, n) {
  total <- 0
  for (j in seq_len(nrow(data))) {
    if (data$f[j] > 0) total <- total + data$f[j] * log(ratio_ab(a, b, data$x[j], n))
  }
  total
}

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
  idx <- seq_len(k + 1)
  for (iter in seq_len(NM_MAX_ITER)) {
    ord <- order(vals, idx)
    pts <- pts[ord]
    vals <- vals[ord]
    size <- 0
    for (i in 2:(k + 1)) size <- max(size, max(abs(pts[[i]] - pts[[1]])))
    if (vals[k + 1] - vals[1] < NM_TOL && size < 1e-8) {
      converged <- TRUE
      break
    }
    centroid <- numeric(k)
    for (j in seq_len(k)) {
      c0 <- 0
      for (i in seq_len(k)) c0 <- c0 + pts[[i]][j]
      centroid[j] <- c0 / k
    }
    refl <- centroid + (centroid - pts[[k + 1]])
    fr <- f(refl)
    if (fr < vals[1]) {
      exp_pt <- centroid + 2 * (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, idx)[1]
  list(x = pts[[best]], value = vals[best], converged = converged)
}

# Maximum likelihood on the log scale, restarted from its own optimum until it stops moving.
fit_bgbb <- function(data, n, start = c(0, 0, 0, 0)) {
  neg <- function(z) if (max(z) > LOG_PARAMETER_MAX) PENALTY else -loglik(exp(z), data, n)
  z <- start
  value <- neg(z)
  for (r in 1:20) {
    res <- nelder_mead(neg, z, 0.5)
    if (!res$converged) fail("the BG/BB fit did not converge within NM_MAX_ITER iterations")
    moved <- value - res$value
    z <- res$x
    value <- res$value
    if (moved < 1e-9) return(list(par = exp(z), ll = -value))
  }
  fail("the BG/BB fit kept moving after 20 restarts")
}

fit_bb <- function(data, n) {
  neg <- function(z) if (max(z) > LOG_PARAMETER_MAX) PENALTY else -bb_loglik(exp(z[1]), exp(z[2]), data, n)
  z <- c(0, 0)
  value <- neg(z)
  for (r in 1:20) {
    res <- nelder_mead(neg, z, 0.5)
    if (!res$converged) fail("the BB fit did not converge within NM_MAX_ITER iterations")
    moved <- value - res$value
    z <- res$x
    value <- res$value
    if (moved < 1e-9) return(list(a = exp(z[1]), b = exp(z[2]), ll = -value))
  }
  fail("the BB fit kept moving after 20 restarts")
}

# Equation (13): expected purchases over the next `horizon` opportunities, as a sum of per-opportunity terms.
expected_next <- function(par, x, tx, n, horizon) {
  head <- ratio_ab(par[1], par[2], x + 1, n + 1) / likelihood(par, x, tx, n)
  total <- 0
  for (k in seq_len(horizon)) total <- total + head * alive_through(par[3], par[4], n + k)
  total
}

# Equation (13) as printed, with gamma functions: an independent route to the same number.
expected_next_closed <- function(par, x, tx, n, horizon) {
  g <- par[3]
  d <- par[4]
  head <- ratio_ab(par[1], par[2], x + 1, n + 1) / likelihood(par, x, tx, n)
  lead <- d / (g - 1) * exp(lgamma(g + d) - lgamma(1 + d))
  bracket <- exp(lgamma(1 + d + n) - lgamma(g + d + n)) - exp(lgamma(1 + d + n + horizon) - lgamma(g + d + n + horizon))
  head * lead * bracket
}

# Equation (11): probability of being alive at opportunity n + 1.
p_alive <- function(par, x, tx, n) {
  ratio_ab(par[1], par[2], x, n) * alive_through(par[3], par[4], n + 1) / likelihood(par, x, tx, n)
}

# Equation (12): probability of xs purchases over the next `horizon` opportunities.
cond_pmf <- function(par, x, tx, n, horizon, xs) {
  lik <- likelihood(par, x, tx, n)
  a2 <- choose(horizon, xs) * ratio_ab(par[1], par[2], x + xs, n + horizon) * alive_through(par[3], par[4], n + horizon)
  if (horizon - 1 >= xs) for (i in xs:(horizon - 1)) {
    a2 <- a2 + choose(i, xs) * ratio_ab(par[1], par[2], x + xs, n + i) * dies_at(par[3], par[4], n + i)
  }
  res <- a2 / lik
  if (xs == 0) res <- res + (1 - ratio_ab(par[1], par[2], x, n) * alive_through(par[3], par[4], n) / lik)
  res
}

# Discounted expected residual transactions: sum over k >= 1 of E[Y_(n+k)] / (1 + rate)^k.
dert <- function(par, x, tx, n, rate, terms = TERMS) {
  g <- par[3]
  d <- par[4]
  head <- ratio_ab(par[1], par[2], x + 1, n + 1) / likelihood(par, x, tx, n)
  term <- alive_through(g, d, n)
  total <- 0
  for (k in seq_len(terms)) {
    term <- term * ((d + n + k - 1) / (g + d + n + k - 1) / (1 + rate))
    total <- total + term
  }
  head * total
}

# Equation (14) of the paper, with 2F1 evaluated by the term recursion of the authors' note.
dert_hypergeometric <- function(par, x, tx, n, rate) {
  g <- par[3]
  d <- par[4]
  z <- 1 / (1 + rate)
  u <- 1
  f21 <- 1
  for (j in seq_len(TERMS)) {
    u <- u * ((1 + j - 1) * (d + n + 1 + j - 1) / ((g + d + n + 1 + j - 1) * j) * z)
    f21 <- f21 + u
  }
  ratio_ab(par[1], par[2], x + 1, n + 1) * alive_through(g, d, n + 1) / (1 + rate) * f21 / likelihood(par, x, tx, n)
}

# Equation (7): probability of x purchases in the first n opportunities.
pmf <- function(par, n, x) {
  res <- choose(n, x) * ratio_ab(par[1], par[2], x, n) * alive_through(par[3], par[4], n)
  if (n - 1 >= x) for (i in x:(n - 1)) res <- res + choose(i, x) * ratio_ab(par[1], par[2], x, i) * dies_at(par[3], par[4], i)
  res
}

# Equation (8).
mean_closed <- function(par, n) {
  a <- par[1]; b <- par[2]; g <- par[3]; d <- par[4]
  a / (a + b) * d / (g - 1) *
    (1 - exp(lgamma(g + d) - lgamma(g + d + n) + lgamma(1 + d + n) - lgamma(1 + d)))
}

# P(x, t_x) = L(x, t_x) times the number of purchase strings with that recency and frequency.
pattern_probabilities <- function(par, n) {
  pt <- patterns(n)
  vapply(seq_len(nrow(pt)), function(j) {
    strings <- if (pt$x[j] == 0) 1 else choose(pt$tx[j] - 1, pt$x[j] - 1)
    strings * likelihood(par, pt$x[j], pt$tx[j], n)
  }, numeric(1))
}

# Identities the reported values rest on; the programs stop if one fails.
check_invariants <- function() {
  par <- c(TRUE_ALPHA, TRUE_BETA, TRUE_GAMMA, TRUE_DELTA)
  n <- CALIBRATION
  probs <- pattern_probabilities(par, n)
  pt <- patterns(n)
  by_x <- numeric(n + 1)
  for (j in seq_along(probs)) by_x[pt$x[j] + 1] <- by_x[pt$x[j] + 1] + probs[j]
  checks <- list(
    list("pattern probabilities sum to one", plain_sum(probs), 1),
    list("pmf sums to one", plain_sum(vapply(0:n, function(x) pmf(par, n, x), numeric(1))), 1),
    list("pmf mean equals equation (8)", plain_sum(vapply(0:n, function(x) x * pmf(par, n, x), numeric(1))), mean_closed(par, n)),
    list("conditional pmf sums to one", plain_sum(vapply(0:HOLDOUT, function(k) cond_pmf(par, 2, 4, n, HOLDOUT, k), numeric(1))), 1),
    list("conditional pmf mean equals equation (13)",
         plain_sum(vapply(0:HOLDOUT, function(k) k * cond_pmf(par, 2, 4, n, HOLDOUT, k), numeric(1))),
         expected_next(par, 2, 4, n, HOLDOUT)),
    list("per-season terms equal the gamma form of equation (13)",
         expected_next(par, 3, 5, n, HOLDOUT), expected_next_closed(par, 3, 5, n, HOLDOUT)),
    list("DERT series equals equation (14)", dert(par, 4, 6, n, DISCOUNT), dert_hypergeometric(par, 4, 6, n, DISCOUNT)),
    list("DERT discounts the first future season once",
         dert(par, 1, 1, n, DISCOUNT, terms = 1L), expected_next(par, 1, 1, n, 1) / (1 + DISCOUNT))
  )
  for (x in 0:n) checks[[length(checks) + 1]] <- list(paste0("patterns with x=", x, " add up to the pmf"), by_x[x + 1], pmf(par, n, x))
  for (check in checks) {
    if (abs(check[[2]] - check[[3]]) > 1e-9 * max(1, abs(check[[3]]))) fail(paste0("invariant failed, ", check[[1]]))
  }
  out("metric.check.invariants=", length(checks), " passed")
}

PUBLISHED <- list(
  # Fader and Hardie (2011), note on the Excel implementation, p. 4: LL at the starting values 1, 1, 1, 1.
  start = "loglik=-37232.0",
  # Fader, Hardie and Shang (2010), Table 4.
  bgbb = "alpha=1.204 beta=0.750 gamma=0.657 delta=2.783 loglik=-33225.6",
  bb = "alpha=0.487 beta=0.826 loglik=-35516.1",
  # Table 5, row by row (x = 0, then x = 1 with t_x = 1..6, ..., x = 6), two decimals.
  table5 = "0.07 | 0.09 0.31 0.59 0.84 1.02 1.15 | 0.12 0.54 1.06 1.44 1.67 | 0.22 1.03 1.80 2.19 | 0.58 2.03 2.71 | 1.81 3.23 | 3.75",
  # Table 6, P(alive in 2002), same layout.
  table6 = "0.11 | 0.07 0.25 0.48 0.68 0.83 0.93 | 0.07 0.30 0.59 0.80 0.93 | 0.10 0.44 0.77 0.93 | 0.20 0.70 0.93 | 0.52 0.93 | 0.93"
)

table_by_row <- function(n, value) {
  rows <- fmt(value(0, 0), 2)
  for (x in 1:n) rows <- c(rows, paste(vapply(x:n, function(tx) fmt(value(x, tx), 2), ""), collapse = " "))
  paste(rows, collapse = " | ")
}

# Reproduce what the sources print before trusting the code on new data; stop if one is missed.
published_checks <- function() {
  n <- 6L
  donors <- plain_sum(DONORS$f)
  repeat_donations <- plain_sum(DONORS$x * DONORS$f)
  if (donors != 11104 || repeat_donations != 24615) fail("Table 2 transcription")
  out("metric.check.table2 patterns=", nrow(DONORS), " donors=", fmt(donors, 0), " repeat_donations=", fmt(repeat_donations, 0))
  line <- paste0("loglik=", fmt(loglik(c(1, 1, 1, 1), DONORS, n), 1))
  if (line != PUBLISHED$start) fail(paste0("Excel note starting log-likelihood not reproduced: ", line))
  out("metric.check.fader_hardie_2011_start ", line)
  fit <- fit_bgbb(DONORS, n)
  par <- fit$par
  line <- paste0("alpha=", fmt(par[1], 3), " beta=", fmt(par[2], 3), " gamma=", fmt(par[3], 3), " delta=", fmt(par[4], 3),
                 " loglik=", fmt(fit$ll, 1))
  if (line != PUBLISHED$bgbb) fail(paste0("Fader, Hardie and Shang (2010) Table 4 BG/BB not reproduced: ", line))
  out("metric.check.fhs_2010_table4_bgbb ", line)
  other <- fit_bgbb(DONORS, n, start = rep(log(0.01), 4))
  if (abs(other$ll - fit$ll) > 1e-6 || max(abs(other$par / par - 1)) > 1e-4) fail("the BG/BB fit depends on its starting values")
  out("metric.check.fhs_2010_table4_bgbb_from_0.01 same_optimum=yes")
  bb <- fit_bb(DONORS, n)
  line <- paste0("alpha=", fmt(bb$a, 3), " beta=", fmt(bb$b, 3), " loglik=", fmt(bb$ll, 1))
  if (line != PUBLISHED$bb) fail(paste0("Fader, Hardie and Shang (2010) Table 4 BB not reproduced: ", line))
  out("metric.check.fhs_2010_table4_bb ", line)
  tables <- list(table5 = function(x, tx) expected_next(par, x, tx, n, 5L),
                 table6 = function(x, tx) p_alive(par, x, tx, n))
  for (name in names(tables)) {
    line <- table_by_row(n, tables[[name]])
    if (line != PUBLISHED[[name]]) fail(paste0("Fader, Hardie and Shang (2010) ", name, " not reproduced: ", line))
    out("metric.check.fhs_2010_", name, "=", line)
  }
  out("metric.check.fhs_2010_zero_repeat_donors_2002_2006=", fmt(3464 * expected_next(par, 0, 0, n, 5L), 1))
  out("metric.check.fhs_2010_prior_means E(P)=", fmt(par[1] / (par[1] + par[2]), 2),
      " E(Theta)=", fmt(par[3] / (par[3] + par[4]), 2))
}

# Exact draws of purchase strings by the Polya-urn form of the model: no beta variate is needed.
# At season t, an alive customer dies with probability g / (g + d + t - 1); if still alive he buys with
# probability (a + purchases so far) / (a + b + seasons alive so far).
simulate_cohort <- function(par, count, seasons) {
  out_m <- matrix(0L, count, seasons)
  for (c0 in seq_len(count)) {
    alive <- TRUE
    bought <- 0
    for (t in seq_len(seasons)) {
      u_die <- uniform()
      u_buy <- uniform()
      if (alive && u_die < par[3] / (par[3] + par[4] + t - 1)) alive <- FALSE
      if (alive && u_buy < (par[1] + bought) / (par[1] + par[2] + t - 1)) {
        out_m[c0, t] <- 1L
        bought <- bought + 1
      }
    }
  }
  out_m
}

recency <- function(ys, n) {
  hits <- which(ys[1:n] == 1L)
  if (length(hits)) max(hits) else 0L
}

# Recency/frequency counts over the first n opportunities, in pattern order.
summarize <- function(strings, n) {
  pt <- patterns(n)
  counts <- integer(nrow(pt))
  for (c0 in seq_len(nrow(strings))) {
    x <- sum(strings[c0, 1:n])
    tx <- recency(strings[c0, ], n)
    j <- which(pt$x == x & pt$tx == tx)
    counts[j] <- counts[j] + 1L
  }
  data.frame(x = pt$x, tx = pt$tx, f = counts)
}

# Multinomial draw of `count` customers over the patterns, by inversion of the cumulative probabilities.
draw_patterns <- function(probs, count) {
  cum <- numeric(length(probs))
  total <- 0
  for (j in seq_along(probs)) {
    total <- total + probs[j]
    cum[j] <- total
  }
  counts <- integer(length(probs))
  last <- length(cum)
  for (i in seq_len(count)) {
    u <- uniform() * total
    j <- 1L
    while (j < last && u >= cum[j]) j <- j + 1L
    counts[j] <- counts[j] + 1L
  }
  counts
}

naive_rate <- function(x, n) x / n

# Inactivity rule: a customer without a purchase in the last two seasons is written off.
rule_rate <- function(x, tx, n) if (n - tx >= 2) 0 else x / n

# Value of a constant expected purchase rate from season n + 1 on, first season discounted once.
perpetuity <- function(rate_per_season, rate = DISCOUNT) MARGIN * rate_per_season / rate

# Damped rhythm: x/n purchases per season, times decay^k at future season k.
damped_purchases <- function(rate_per_season, decay, horizon) {
  total <- 0
  factor <- 1
  for (k in seq_len(horizon)) {
    factor <- factor * decay
    total <- total + factor
  }
  rate_per_season * total
}

# Value of the damped rhythm: m (x/n) sum_k decay^k / (1 + rate)^k = m (x/n) decay / (1 + rate - decay).
damped_value <- function(rate_per_season, decay, rate = DISCOUNT) MARGIN * rate_per_season * decay / (1 + rate - decay)

# P(chi2_df > stat) = 1 - P(df/2, stat/2), the regularized lower gamma by its power series (as in Python, not pchisq).
chi2_upper_tail <- function(stat, df) {
  a <- df / 2
  x <- stat / 2
  term <- 1 / a
  total <- term
  for (k in 1:999) {
    term <- term * (x / (a + k))
    total <- total + term
    if (term < 1e-17 * total) break
  }
  1 - exp(a * log(x) - x - lgamma(a)) * total
}

# Smallest k with P(X <= k | history) >= 1/2 over the next `horizon` opportunities.
predictive_median <- function(par, x, tx, n, horizon) {
  cumulative <- 0
  for (k in 0:horizon) {
    cumulative <- cumulative + cond_pmf(par, x, tx, n, horizon, k)
    if (cumulative >= 0.5) return(k)
  }
  horizon
}

true <- c(TRUE_ALPHA, TRUE_BETA, TRUE_GAMMA, TRUE_DELTA)
n <- CALIBRATION
out("metric.design.customers=", CUSTOMERS, " calibration_seasons=", n, " holdout_seasons=", HOLDOUT)
out("metric.design.seed=", LCG_SEED, " bootstrap=", BOOTSTRAP, " replications=", REPLICATIONS)
out("metric.design.true_alpha=", fmt(TRUE_ALPHA, 2), " true_beta=", fmt(TRUE_BETA, 2), " true_gamma=", fmt(TRUE_GAMMA, 2),
    " true_delta=", fmt(TRUE_DELTA, 2), " true_mean_p=", fmt(TRUE_ALPHA / (TRUE_ALPHA + TRUE_BETA), 4),
    " true_mean_theta=", fmt(TRUE_GAMMA / (TRUE_GAMMA + TRUE_DELTA), 4))
out("metric.design.margin_per_pass=", fmt(MARGIN, 2), " discount_per_season=", fmt(DISCOUNT, 4), " series_terms=", TERMS)
published_checks()
check_invariants()

strings <- simulate_cohort(true, CUSTOMERS, n + HOLDOUT)
data <- summarize(strings, n)
holdout <- vapply(seq_len(CUSTOMERS), function(c0) sum(strings[c0, (n + 1):(n + HOLDOUT)]), integer(1))
out("metric.data.patterns=", paste0(data$x, "/", data$tx, ":", data$f, collapse = " "))
buyers <- vapply(seq_len(n + HOLDOUT), function(t) sum(strings[, t]), integer(1))
out("metric.data.buyers_by_season=", paste0(seq_len(n + HOLDOUT), ":", buyers, collapse = " "))
zero <- sum(data$f[data$x == 0])
lapsed <- sum(data$f[n - data$tx >= 2])
out("metric.data.calibration_repeat_purchases=", sum(data$x * data$f), " holdout_purchases=", sum(holdout),
    " never_returned=", zero, " lapsed_two_seasons_or_more=", lapsed)
out("metric.data.share never_returned=", fmt(100 * zero / CUSTOMERS, 1), "%",
    " lapsed_two_seasons_or_more=", fmt(100 * lapsed / CUSTOMERS, 1), "%")
# The damped shortcut lets the past rhythm decay at the mean rate at which the number of buyers holds up, seasons 1 to n.
decay <- (buyers[n] / buyers[1])^(1 / (n - 1))
out("metric.data.buyer_count_ratio_per_season seasons_1_to_", n, "=", fmt(decay, 3))

fit <- fit_bgbb(data, n)
par <- fit$par
other <- fit_bgbb(data, n, start = rep(log(0.01), 4))
if (abs(other$ll - fit$ll) > 1e-6 || max(abs(other$par / par - 1)) > 1e-4) fail("the synthetic BG/BB fit depends on its starting values")
out("metric.check.synthetic_fit_from_0.01 same_optimum=yes")
bb <- fit_bb(data, n)
a_bb <- bb$a
b_bb <- bb$b
out("metric.fit.bgbb alpha=", fmt(par[1], 3), " beta=", fmt(par[2], 3), " gamma=", fmt(par[3], 3), " delta=", fmt(par[4], 3),
    " loglik=", fmt(fit$ll, 1))
out("metric.fit.bgbb mean_p=", fmt(par[1] / (par[1] + par[2]), 4), " mean_theta=", fmt(par[3] / (par[3] + par[4]), 4))
out("metric.fit.bb alpha=", fmt(a_bb, 3), " beta=", fmt(b_bb, 3), " loglik=", fmt(bb$ll, 1),
    " likelihood_ratio=", fmt(2 * (fit$ll - bb$ll), 1))
expected_counts <- vapply(0:n, function(x) CUSTOMERS * pmf(par, n, x), numeric(1))
actual_counts <- vapply(0:n, function(x) sum(data$f[data$x == x]), integer(1))
chi2 <- plain_sum((actual_counts - expected_counts)^2 / expected_counts)
out("metric.fit.frequency_actual_vs_bgbb=",
    paste0(0:n, ":", actual_counts, "/", vapply(expected_counts, function(v) fmt(v, 1), ""), collapse = " "),
    " chi2=", fmt(chi2, 2))
# Fit on the 22 recency/frequency cells the model is estimated on: 22 - 1 - 4 parameters = 17 degrees of freedom.
cell_expected <- CUSTOMERS * pattern_probabilities(par, n)
chi2_cells <- plain_sum((data$f - cell_expected)^2 / cell_expected)
df_cells <- nrow(data) - 1L - 4L
out("metric.fit.chi2_22_cells chi2=", fmt(chi2_cells, 2), " df=", df_cells, " p=", fmt(chi2_upper_tail(chi2_cells, df_cells), 2),
    " min_expected=", fmt(min(cell_expected), 1))

# Holdout forecasts per customer: four methods and the expectation under the true parameters.
methods <- c("naive", "rule", "damped", "bb", "bgbb", "truth")
per_pattern <- t(vapply(seq_len(nrow(data)), function(j) {
  x <- data$x[j]
  tx <- data$tx[j]
  c(naive = HOLDOUT * naive_rate(x, n), rule = HOLDOUT * rule_rate(x, tx, n),
    damped = damped_purchases(naive_rate(x, n), decay, HOLDOUT),
    bb = HOLDOUT * (a_bb + x) / (a_bb + b_bb + n),
    bgbb = expected_next(par, x, tx, n, HOLDOUT), truth = expected_next(true, x, tx, n, HOLDOUT))
}, numeric(6)))
key_x <- vapply(seq_len(CUSTOMERS), function(c0) sum(strings[c0, 1:n]), integer(1))
key_tx <- vapply(seq_len(CUSTOMERS), function(c0) recency(strings[c0, ], n), integer(1))
key_row <- vapply(seq_len(CUSTOMERS), function(c0) which(data$x == key_x[c0] & data$tx == key_tx[c0]), integer(1))
totals <- vapply(methods, function(m) plain_sum(per_pattern[key_row, m]), numeric(1))
out("metric.holdout.total actual=", sum(holdout), " ",
    paste0(methods, "=", vapply(totals, function(v) fmt(v, 1), ""), collapse = " "))
out("metric.holdout.total_gap_vs_actual ",
    paste0(methods, "=", vapply(totals, function(v) fmt(100 * (v / sum(holdout) - 1), 1), ""), "%", collapse = " "))
# Chance alone: standard deviation of the realized total given the histories, under the true parameters.
variance <- 0
for (c0 in seq_len(CUSTOMERS)) {
  mean_k <- 0
  second <- 0
  for (xs in 0:HOLDOUT) {
    pr <- cond_pmf(true, key_x[c0], key_tx[c0], n, HOLDOUT, xs)
    mean_k <- mean_k + xs * pr
    second <- second + xs * xs * pr
  }
  variance <- variance + (second - mean_k * mean_k)
}
out("metric.holdout.chance_sd_of_total truth_sd=", fmt(sqrt(variance), 1),
    " share_of_truth=", fmt(100 * sqrt(variance) / totals[["truth"]], 1), "%")
for (label in c("frequency", "recency")) {
  group <- if (label == "frequency") key_x else key_tx
  for (v in 0:n) {
    members <- which(group == v)
    if (!length(members)) next
    m_actual <- sum(holdout[members]) / length(members)
    out("metric.holdout.by_", label, "=", v, " customers=", length(members), " actual=", fmt(m_actual, 2), " ",
        paste0(methods, "=", vapply(methods, function(m) fmt(plain_sum(per_pattern[key_row[members], m]) / length(members), 2), ""),
               collapse = " "))
  }
}
mae <- vapply(methods, function(m) plain_sum(abs(per_pattern[key_row, m] - holdout)) / CUSTOMERS, numeric(1))
out("metric.holdout.mean_abs_error_per_customer ", paste0(methods, "=", vapply(mae, function(v) fmt(v, 3), ""), collapse = " "))
medians <- vapply(seq_len(nrow(data)), function(j) predictive_median(par, data$x[j], data$tx[j], n, HOLDOUT), numeric(1))
mae_median <- plain_sum(abs(medians[key_row] - holdout)) / CUSTOMERS
out("metric.holdout.mean_abs_error_bgbb_predictive_median=", fmt(mae_median, 3))
active_pred <- plain_sum(vapply(seq_len(CUSTOMERS), function(c0) 1 - cond_pmf(par, key_x[c0], key_tx[c0], n, HOLDOUT, 0), numeric(1)))
active_true <- plain_sum(vapply(seq_len(CUSTOMERS), function(c0) 1 - cond_pmf(true, key_x[c0], key_tx[c0], n, HOLDOUT, 0), numeric(1)))
out("metric.holdout.active_customers actual=", sum(holdout > 0), " bgbb=", fmt(active_pred, 1), " truth=", fmt(active_true, 1))

for (profile in list(c(6, 6), c(5, 5), c(4, 6), c(3, 3), c(1, 6), c(0, 0))) {
  x <- profile[1]
  tx <- profile[2]
  j <- which(data$x == x & data$tx == tx)
  out("metric.profile.x=", x, " tx=", tx, " customers=", data$f[j],
      " naive=", fmt(per_pattern[j, "naive"], 2), " rule=", fmt(per_pattern[j, "rule"], 2), " bb=", fmt(per_pattern[j, "bb"], 2),
      " bgbb=", fmt(per_pattern[j, "bgbb"], 2), " truth=", fmt(per_pattern[j, "truth"], 2),
      " p_alive_bgbb=", fmt(p_alive(par, x, tx, n), 2), " p_alive_truth=", fmt(p_alive(true, x, tx, n), 2),
      " value_bgbb=", fmt(MARGIN * dert(par, x, tx, n, DISCOUNT), 2),
      " value_truth=", fmt(MARGIN * dert(true, x, tx, n, DISCOUNT), 2),
      " value_naive=", fmt(perpetuity(naive_rate(x, n)), 2))
}

# Values of the cohort: residual lifetime value, first future season discounted once.
cohort_value <- function(fn) plain_sum(vapply(seq_len(nrow(data)), function(j) data$f[j] * fn(data$x[j], data$tx[j]), numeric(1)))
values <- c(
  naive = cohort_value(function(x, tx) perpetuity(naive_rate(x, n))),
  rule = cohort_value(function(x, tx) perpetuity(rule_rate(x, tx, n))),
  damped = cohort_value(function(x, tx) damped_value(naive_rate(x, n), decay)),
  bb = cohort_value(function(x, tx) perpetuity((a_bb + x) / (a_bb + b_bb + n))),
  bgbb = cohort_value(function(x, tx) MARGIN * dert(par, x, tx, n, DISCOUNT)),
  truth = cohort_value(function(x, tx) MARGIN * dert(true, x, tx, n, DISCOUNT))
)
out("metric.value.cohort ", paste0(methods, "=", vapply(values, function(v) fmt(v, 0), ""), collapse = " "))
out("metric.value.cohort_gap_vs_truth ",
    paste0(methods, "=", vapply(values, function(v) fmt(100 * (v / values[["truth"]] - 1), 1), ""), "%", collapse = " "))
out("metric.value.per_customer ", paste0(methods, "=", vapply(values, function(v) fmt(v / CUSTOMERS, 3), ""), collapse = " "))
lapsed_rows <- which(n - data$tx >= 2)
lapsed_value <- plain_sum(vapply(lapsed_rows, function(j) data$f[j] * MARGIN * dert(par, data$x[j], data$tx[j], n, DISCOUNT), numeric(1)))
lapsed_customers <- which(n - key_tx >= 2)
lapsed_holdout <- sum(holdout[lapsed_customers])
lapsed_active <- sum(holdout[lapsed_customers] > 0)
lapsed_expected <- plain_sum(per_pattern[key_row[lapsed_customers], "bgbb"])
out("metric.value.written_off_by_rule customers=", lapsed, " bgbb_value=", fmt(lapsed_value, 0),
    " share_of_bgbb_value=", fmt(100 * lapsed_value / values[["bgbb"]], 1), "%",
    " holdout_purchases_actual=", lapsed_holdout, " holdout_purchases_bgbb=", fmt(lapsed_expected, 1),
    " customers_back_in_holdout=", lapsed_active)
zero_value <- plain_sum(vapply(which(data$x == 0), function(j) data$f[j] * MARGIN * dert(par, data$x[j], data$tx[j], n, DISCOUNT), numeric(1)))
out("metric.value.never_returned customers=", zero, " value_each=", fmt(MARGIN * dert(par, 0, 0, n, DISCOUNT), 2),
    " value_total=", fmt(zero_value, 0), " share_of_bgbb_value=", fmt(100 * zero_value / values[["bgbb"]], 1), "%")

cohort_within <- function(seasons) cohort_value(function(x, tx) MARGIN * dert(par, x, tx, n, DISCOUNT, terms = seasons))
out("metric.horizon.share_of_bgbb_value_beyond_season=",
    paste0(n + c(5, 10, 20), ":", vapply(c(5L, 10L, 20L), function(h) fmt(100 * (1 - cohort_within(h) / values[["bgbb"]]), 1), ""),
           "%", collapse = " "))
for (rate in c(0.05, 0.10, 0.15)) {
  v_b <- cohort_value(function(x, tx) MARGIN * dert(par, x, tx, n, rate))
  v_t <- cohort_value(function(x, tx) MARGIN * dert(true, x, tx, n, rate))
  v_n <- cohort_value(function(x, tx) perpetuity(naive_rate(x, n), rate))
  v_r <- cohort_value(function(x, tx) perpetuity(rule_rate(x, tx, n), rate))
  v_d <- cohort_value(function(x, tx) damped_value(naive_rate(x, n), decay, rate))
  out("metric.sensitivity.discount=", fmt(100 * rate, 0), "% cohort_bgbb=", fmt(v_b, 0), " cohort_truth=", fmt(v_t, 0),
      " cohort_naive=", fmt(v_n, 0), " cohort_rule=", fmt(v_r, 0), " cohort_damped=", fmt(v_d, 0),
      " best_customer_bgbb=", fmt(MARGIN * dert(par, 6, 6, n, rate), 2))
}

probs <- pattern_probabilities(par, n)
pt <- patterns(n)
boot_holdout <- numeric(BOOTSTRAP)
boot_value <- numeric(BOOTSTRAP)
boot_best <- numeric(BOOTSTRAP)
boot_mean_p <- numeric(BOOTSTRAP)
boot_mean_theta <- numeric(BOOTSTRAP)
for (b in seq_len(BOOTSTRAP)) {
  counts <- draw_patterns(probs, CUSTOMERS)
  b_data <- data.frame(x = pt$x, tx = pt$tx, f = counts)
  b_par <- fit_bgbb(b_data, n, start = log(par))$par
  boot_holdout[b] <- cohort_value(function(x, tx) expected_next(b_par, x, tx, n, HOLDOUT))
  boot_value[b] <- cohort_value(function(x, tx) MARGIN * dert(b_par, x, tx, n, DISCOUNT))
  boot_best[b] <- MARGIN * dert(b_par, 6, 6, n, DISCOUNT)
  boot_mean_p[b] <- b_par[1] / (b_par[1] + b_par[2])
  boot_mean_theta[b] <- b_par[3] / (b_par[3] + b_par[4])
}
boot_holdout <- sort(boot_holdout)
boot_value <- sort(boot_value)
boot_best <- sort(boot_best)
boot_mean_p <- sort(boot_mean_p)
boot_mean_theta <- sort(boot_mean_theta)
lo <- (BOOTSTRAP * 25L) %/% 1000L
hi <- (BOOTSTRAP * 975L) %/% 1000L
out("metric.uncertainty.parametric_bootstrap draws=", BOOTSTRAP, " sorted_ranks=", lo, "..", hi,
    " holdout_total_95=", fmt(boot_holdout[lo], 1), "..", fmt(boot_holdout[hi], 1),
    " cohort_value_95=", fmt(boot_value[lo], 0), "..", fmt(boot_value[hi], 0),
    " best_customer_value_95=", fmt(boot_best[lo], 2), "..", fmt(boot_best[hi], 2))
out("metric.uncertainty.parametric_bootstrap mean_p_95=", fmt(boot_mean_p[lo], 4), "..", fmt(boot_mean_p[hi], 4),
    " mean_theta_95=", fmt(boot_mean_theta[lo], 4), "..", fmt(boot_mean_theta[hi], 4))

true_probs <- pattern_probabilities(true, n)
# The damped shortcut needs season-by-season buyers, which a recency/frequency draw does not give: it is
# judged on the detailed cohort above, not in the replications.
mc_methods <- c("naive", "rule", "bb", "bgbb", "truth")
rep <- matrix(0, REPLICATIONS, length(mc_methods), dimnames = list(NULL, mc_methods))
naive_above <- 0L
rule_below <- 0L
for (r in seq_len(REPLICATIONS)) {
  counts <- draw_patterns(true_probs, CUSTOMERS)
  r_data <- data.frame(x = pt$x, tx = pt$tx, f = counts)
  r_par <- fit_bgbb(r_data, n)$par
  rbb <- fit_bb(r_data, n)
  over <- function(fn) plain_sum(vapply(seq_len(nrow(r_data)), function(j) r_data$f[j] * fn(r_data$x[j], r_data$tx[j]), numeric(1)))
  rep[r, "naive"] <- over(function(x, tx) HOLDOUT * naive_rate(x, n))
  rep[r, "rule"] <- over(function(x, tx) HOLDOUT * rule_rate(x, tx, n))
  rep[r, "bb"] <- over(function(x, tx) HOLDOUT * (rbb$a + x) / (rbb$a + rbb$b + n))
  rep[r, "bgbb"] <- over(function(x, tx) expected_next(r_par, x, tx, n, HOLDOUT))
  rep[r, "truth"] <- over(function(x, tx) expected_next(true, x, tx, n, HOLDOUT))
  if (rep[r, "naive"] > rep[r, "truth"]) naive_above <- naive_above + 1L
  if (rep[r, "rule"] < rep[r, "truth"]) rule_below <- rule_below + 1L
}
for (m in c("naive", "rule", "bb", "bgbb")) {
  gaps <- rep[, m] / rep[, "truth"] - 1
  mean_gap <- plain_sum(gaps) / REPLICATIONS
  sd_gap <- sqrt(plain_sum((gaps - mean_gap) * (gaps - mean_gap)) / (REPLICATIONS - 1))
  out("metric.monte_carlo.holdout_total_gap_vs_truth ", m, " mean=", fmt(100 * mean_gap, 1), "%",
      " sd=", fmt(100 * sd_gap, 1), "% mcse=", fmt(100 * sd_gap / sqrt(REPLICATIONS), 1), "%",
      " min=", fmt(100 * min(gaps), 1), "% max=", fmt(100 * max(gaps), 1), "%")
}
out("metric.monte_carlo.naive_above_truth=", naive_above, " rule_below_truth=", rule_below, " of ", REPLICATIONS)