ripr

ripr computes the reverse information projection (RIPr) P* of an alternative distribution Q onto a (possibly non-convex, union-structured) null hypothesis, together with the duality-gap certificate for the rescaled e-variable (Q / P*) / (1 + gap).

How strong that certificate is depends on how the gap was obtained, and every certificate says which case it is in via gap_certified. For a multinomial family on simplex-shaped null faces with an exact engine, pass bnb = bnb_control() and the gap comes from a Bernstein branch-and-bound upper bound on sup_theta G: the rescaled e-variable’s expectation under every null distribution is then at most 1 as a proven inequality. Otherwise the gap comes from a multi-start BFGS search, which returns a lower bound on the face maximum — so a missed global optimum leaves the e-variable under-corrected, and the guarantee is evidence rather than proof. With a Monte Carlo engine there is a second, separate error: the gap is estimated, so the bound is stochastic and holds at level conf (default 0.95), at the worst-case theta* the certifying sample selects.

The projection is found by Frank–Wolfe and EM over a mixture of atoms on the null faces. It requires the specification of a sampling family (the model p_theta), the null geometry (a union of faces), and a surrogate alternative Q. The alternative is an outcome_distribution — a law over the sample space — which you can supply directly or build from a mixing measure over the parameter space (e.g. point_mixing, finite_mixing), marginalised through the family. These components are templated as S7 base classes, so new sample spaces drop in without touching the algorithmic core. The algorithms require the ability to compute KL divergences and gradients, which may be implemented either via exact numerical methods or by Monte Carlo sampling. These two routes are encapsulated in the engine class.

# install from github
renv::install("fleverest/ripr")

Example 1: a one-sided binomial null

Let X ~ Binomial(10, p), written as a 2-category multinomial with theta = (p, 1 - p). Consider testing H_0: p <= 1/2 against H_1: p > 1/2, adopting the simple alternative Q = Binomial(10, 3/4). The RIPr of Q onto the convex null is a point mass at p = 1/2. We demonstrate this using both the exact and Monte Carlo engines.

library(ripr)
n <- 10
fam <- multinomial_family(n_trials = n, k = 2)
# `ripr_problem` takes Q as an `outcome_distribution` over the sample space.
# Build one by marginalising a `mixing` measure over the parameter theta (here
# it is just a point mass at theta = (3/4, 1/4)) through the family, giving the
# outcome law Q = Binomial(10, 3/4).
Q <- as_marginal(point_mixing(theta_star = c(0.75, 0.25)), fam)

# H_0 = {p <= 1/2} is the simplex segment from p = 0 to the tie p = 1/2.
face <- polytope_face(vertices = cbind(c(0, 1), c(0.5, 0.5)), face_index = 1)
null <- null_region(faces = list(face))

Exact engine (enumerates all 11 outcomes):

set.seed(1)
prob <- ripr_problem(fam, null, Q)
res <- run_ripr(
  prob,
  # Deliberately initialising somewhere sub-optimal for illustrative purposes
  init_atoms = matrix(c(0.25, 0.75), ncol = 1),
  init_atom_faces = 1L,
  fw_iters = 8,
  em_iters = 5,
  prune_threshold = 1e-6, # drop the numerically-dead atoms for output
  verbose = FALSE
)

cbind(
  p = res$projection@mixing@components[1, ],
  weight = res$projection@mixing@weights
)
##              p weight
## [1,] 0.4999973      1

res$projection is a marginal — the fitted P* as a distribution over outcomes — and its @mixing is the finite_mixing holding the atoms and weights on the null. run_ripr() returns the fitted e-variable directly; call e_value() on new data. Observing (8, 2) yields evidence against H_0, while a tie (5, 5) yields evidence in favour of it:

res$e_variable
## <e_variable>  e(x) = (Q / P*) / (1 + gap)
##   numerator  Q  : point_mixing
##   projection P* : finite_mixing (1 atom)
##   gap           : 0   (correction 1 + gap = 1)
##   gap certified : FALSE  (heuristic oracle; see ?certify)

e_value(res$e_variable, rbind(c(8, 2), c(5, 5)))
## [1] 6.4074351 0.2373047

Certifying the gap deterministically

gap is what makes the rescaling valid, and by default it comes from a search: oracle() maximises G(theta) = E_Q[p_theta / P*] by multi-start BFGS, so it returns a lower bound on the face maximum. For a multinomial family on simplex-shaped faces with an exact engine, the gap can instead be proved. p_theta(x) is the degree-n Bernstein basis function at the count vector x, so G already is a Bernstein form, with coefficients q(x) / P*(x). That basis is non-negative and sums to one, so the largest coefficient bounds G across the whole face for free, and de Casteljau subdivision refines it. certifiable() reports whether a problem qualifies:

certifiable(face, prob$engine, fam)
## [1] TRUE

heuristic <- certify(res$projection, prob, n_draws = 1e5)
certified <- certify(res$projection, prob, bnb = bnb_control())

data.frame(
  method = c(heuristic$gap_method, certified$gap_method),
  theta_star = c(heuristic$oracle_theta[1], certified$oracle_theta[1]),
  gap_used = c(heuristic$gap_used, certified$gap_used),
  certified = c(heuristic$gap_certified, certified$gap_certified)
)
##            method theta_star     gap_used certified
## 1 multistart_bfgs  0.4999974 1.135024e-06     FALSE
## 2   bernstein_bnb  0.5000000 2.712102e-05      TRUE

The certified gap is the larger of the two, and that is the whole point. G is increasing on this null, so its supremum sits at the face’s boundary vertex p = 1/2 — but oracle() searches in a softmax parametrisation that covers only the face’s interior, and the objective is nearly flat there, so BFGS stops just short, at the projection’s own atom. The bound evaluates the vertex exactly.

Small as the difference is (in this example), it is the difference between an e-variable and not one. E_theta[e(X)] must be at most 1 at every theta inthe null:

X <- support(fam)
null_expectation <- function(gap, p) {
  ev <- e_variable(numerator = Q, projection = res$projection, gap = gap)
  sum(exp(log_density(fam, c(p, 1 - p)) + e_value(ev, X, log = TRUE)))
}
# Distance above 1 at the worst-case theta. Positive is a violation.
c(
  heuristic = null_expectation(heuristic$gap_used, 0.5) - 1,
  certified = null_expectation(certified$gap_used, 0.5) - 1
)
##     heuristic     certified 
##  2.598597e-05 -1.332268e-15

bnb_iterations is 0 here: the largest Bernstein coefficient already sits at the vertex attaining the maximum; with a single atom mixture, the likelihood ratio in monotone in p. The bound is exactly tight before any subdivision happens. Harder faces bisect until bound - incumbent <= tol. The bound is valid at every iteration, so tol and max_iter buy tightness, never validity.

To carry the guarantee through a fit rather than re-certifying afterwards, pass certify_bnb — the returned e_variable then records the guarantee level, so a caller holding only that object can still tell a proven rescaling from a heuristic one:

set.seed(1)
res_cert <- run_ripr(
  prob,
  init_atoms = matrix(c(0.25, 0.75), ncol = 1),
  init_atom_faces = 1L,
  fw_iters = 8,
  em_iters = 5,
  prune_threshold = 1e-6,
  verbose = FALSE,
  certify_bnb = bnb_control()
)
res_cert$e_variable
## <e_variable>  e(x) = (Q / P*) / (1 + gap)
##   numerator  Q  : point_mixing
##   projection P* : finite_mixing (1 atom)
##   gap           : 2.712e-05   (correction 1 + gap = 1.00003)
##   gap certified : TRUE

This path needs all three of an exact engine, a family whose density is a Bernstein basis function of theta (the multinomial’s is; a Gaussian’s is not), and bounded simplex faces — halfspace_face has no simplex to subdivide, which is a permanent gap rather than a missing feature. Every other combination degrades to the heuristic path and reports gap_certified = FALSE rather than pretending; require_certified = TRUE turns that fallback into an error naming what blocked it. See ?certify for what the two regimes cover, and bnb_control()’s max_coef for the cost ceiling — coefficients per sub-simplex are choose(n + K - 1, K - 1), so certification caps the usable batch size more tightly than optimisation does.

Monte Carlo engine (1000 draws from Q) reaches the same projection:

set.seed(3)
eng_mc <- mc_engine(fam, Q, n_draws = 1000)
prob_mc <- ripr_problem(fam, null, Q, engine = eng_mc)
res_mc <- run_ripr(
  prob_mc,
  init_atoms = matrix(c(0.25, 0.75), ncol = 1),
  init_atom_faces = 1L,
  fw_iters = 10,
  em_iters = 10,
  gap_tol = 1e-3,
  verbose = FALSE
)

cbind(
  p = res_mc$projection@mixing@components[1, ],
  weight = res_mc$projection@mixing@weights
)
##              p weight
## [1,] 0.4999895      1

set.seed() matters here: an mc_engine freezes its draw set at construction, so a run is reproducible, but certification resamples and is not deterministic across seeds.

The result of a fit

run_ripr() returns six elements, partitioned by where each number came from: projection and e_variable are the deliverables, certificate holds everything measured on a fresh certification sample, history and checkpoints hold everything measured on the fit sample, and converged reports whether the fit gap met gap_tol.

names(res_mc)
## [1] "projection"  "e_variable"  "certificate" "history"     "checkpoints"
## [6] "converged"

The certificate reports the guaranteed e-value growth rate, the standard error of the estimated gap, and the settings that produced it. Certification always runs on draws independent of the fit, and on 10x as many by default – it is a single oracle sweep, not a per-iteration cost, so it is cheap to certify far more thoroughly than you fit:

unlist(res_mc$certificate[c(
  "gap",
  "gap_se",
  "gap_used",
  "growth_rate",
  "n_draws",
  "ess"
)])
##          gap       gap_se     gap_used  growth_rate      n_draws          ess 
## 8.063488e-05 4.393297e-07 8.135751e-05 1.338687e+00 1.000000e+04 1.000000e+04

history is one row per outer iteration, with the inner init/FW/EM steps and the oracle argmax nested as list columns:

res_mc$history[, c("iter", "gap", "gap_se", "support_size", "kl_after_em")]
##   iter          gap       gap_se support_size kl_after_em
## 1    0 7.706625e-05 1.409271e-06            1    1.234616
head(res_mc$history$kl_trace[[1]], 4) # the init + EM steps inside iteration 0
##   step_type n_atoms       kl
## 1      init       1 5.345847
## 2        em       1 1.234616
## 3        em       1 1.234616

checkpoints$final is always present and describes the returned projection after pruning. Its three fields are exactly the init_atoms, init_atom_faces and init_weights arguments needed to resume the fit, so a run can be continued without re-deriving anything (init_weights matters once a projection has more than one atom — without it the atoms carry over but the mixture restarts from uniform weights):

str(res_mc$checkpoints$final, max.level = 1)
## List of 5
##  $ iter         : int NA
##  $ atoms        :List of 1
##  $ weights      : num 1
##  $ atom_face_idx: int 1
##  $ oracle_theta : logi NA

Re-certifying on a larger sample

The Monte Carlo certificate is a stochastic bound, and deliberately a conservative one. certify() resamples fresh draws from Q (it does not reuse the fit draws) and takes the oracle maximum over the null faces; that maximum is biased upward for the true sup_theta G, which is the safe direction, and gap_used inflates it further to a one-sided bound at level conf. gap_used is the only quantity that may be used to rescale an e-variable. Re-certifying the same fitted projection on more draws needs no re-fitting — just raise n_draws:

set.seed(9)
cert <- certify(res_mc$projection, prob_mc, n_draws = 1e5)
c(gap = cert$gap, se = cert$gap_se, gap_used = cert$gap_used, ess = cert$ess)
##          gap           se     gap_used          ess 
## 7.979517e-05 1.377463e-07 8.002174e-05 1.000000e+05

certify() takes the fitted projection (or any finite_mixing of atoms on the null) and sweeps the oracle over the null faces. It resamples the problem’s own engine by default; pass engine = to certify against a different one (say an exact_engine).

estimate = TRUE draws a second independent sample and adds gap_est, an unbiased point estimate of the gap at the selected theta*. It is a diagnostic — useful for comparing fw_variants or deciding whether to keep optimising — and is biased downward relative to sup_theta G. It must never be used to rescale an e-variable.

Example 2: A 4-category plurality null

Let X ~ Multinomial(25, theta). Consider testing H_0: theta_1 <= theta_i for some i!=1 against H_1: theta_1 > theta_i for all i!=1, adopting the simple alternative Q = Multinomial(25, (0.3, 0.27, 0.23, 0.15, 0.05)). The RIPr of Q onto the null is no longer a point mass on the boundary.

n <- 25
k <- 5
q <- c(0.3, 0.27, 0.23, 0.15, 0.05)

# The plurality null can be written as H_0 = H_02 ∪ H_03 ∪ H_04, where
# H_0i = {theta_1 <= theta_i} is a subsimplex with vertices satisfying
# (i) theta_1=theta_i=1/2, (ii) theta_1=theta_i=0, and
# (iii) theta_j = 1 for all other j.
plurality_faces <- function(K) {
  lapply(2:K, function(j) {
    basis <- lapply(setdiff(seq_len(K), 1L), function(k) {
      v <- numeric(K)
      v[k] <- 1
      v
    })
    tie <- numeric(K)
    tie[c(1L, j)] <- 0.5
    polytope_face(
      vertices = do.call(cbind, c(basis, list(tie))),
      face_index = j
    )
  })
}

fam <- multinomial_family(n_trials = n, k = k)
Q <- as_marginal(point_mixing(theta_star = q), fam)
faces <- plurality_faces(k)
null <- null_region(faces = faces)
prob <- ripr_problem(fam, null, Q)

# atoms to start from: one for each face
init <- do.call(
  cbind,
  lapply(faces, function(f) init_point(f, q))
)

res <- run_ripr(
  prob,
  init_atoms = init,
  init_atom_faces = seq_len(k - 1),
  fw_iters = 10,
  em_iters = 10,
  verbose = FALSE
)

heuristic <- certify(res$projection, prob, n_draws = 1e5)
certified <- certify(res$projection, prob, bnb = bnb_control(tol = 1e-9))

data.frame(
  method = c(heuristic$gap_method, certified$gap_method),
  gap_used = c(heuristic$gap_used, certified$gap_used),
  certified = c(heuristic$gap_certified, certified$gap_certified)
)
##            method    gap_used certified
## 1 multistart_bfgs 0.002688946     FALSE
## 2   bernstein_bnb 0.002763226      TRUE

Example 3

TODO: Gaussian plurality via quadrature/monte carlo.