Getting Started with ripr

library(ripr)

Introduction

ripr targets the growth-rate optimal e-variable of safe testing (Grünwald et al. 2024), approximating the reverse information projection (RIPr) of a fixed alternative Q onto a composite null H0 that takes the form of a union of convex parts.

The idea it rests on is simpler than the projection itself. For any non-negative, bounded random variable X, if b is any number at least sup_{theta in Theta_0} E_theta[X], then X / b has null expectation at most one. Specifically, E = X / b is an e-variable.

This means the projection need not be fit exactly: any approximation can be used, so long as the resulting likelihood ratio is corrected against a proven bound to yield a valid e-variable. The quality of the approximation determines only how much e-power survives the correction. That splits a hard, non-convex optimisation problem into an easy one (approximate) and a proof (bound), and this split is what ripr facilitates: fitting is handled by sliding Frank–Wolfe, and certification by Bernstein branch-and-bound.

This vignette walks through both halves on a worked example, then shows what happens when the fit is deliberately bad, and finally contrasts a proven bound against an uncertified search.

A first null: the plurality null

Take a three-category multinomial distribution (n = 20) and the plurality null: candidate 1 does not win outright, H0 = union_j {theta : theta_1 <= theta_j}. This null is a union of convex parts (one region where theta_1 <= theta_2 and another where theta_1 <= theta_3, each part a simplex), so the plurality null is a single null_model holding two simplex_region parts.

K <- 3L
n <- 20L
family <- multinomial_family(n_trials = n, k = K)

# Defines {theta : theta_1 <= theta_j}
plurality_part <- function(j) {
  vertices <- diag(K)
  vertices[1L, ] <- replace(numeric(K), c(1L, j), 0.5)
  simplex_region(vertices = vertices)
}


plurality <- null_model(
  family,
  plurality_part(2) | plurality_part(3)
)

We nominate a discrete mixture-multinomial distribution as the alternative. Here the mixing measure is a single dirac at q = (0.40, 0.35, 0.25), yielding a simple Multinomial(20, q) alternative.

Q <- family(c(0.40, 0.35, 0.25))

Plurality null, with the support of the alternative mixing distribution W1.

The two parts overlap, and the alternative sits in the cap above them, where candidate 1 leads.

Fitting the projection

The fitting procedure is a sequence of steps: e.g. sfw_step() for sliding Frank–Wolfe updates and em_step() for EM (also fw_step for standard Frank–Wolfe variants, and lb_step() for Li–Barron (Li and Barron 1999), though Li–Barron is not generally recommended: see ?oracles). There is no fixed pipeline; which algorithm runs is decided by which verbs are called, and in what order. Here we run 40 steps of Frank–Wolfe, then use ripr_finish() to solve the weights and drop any low-mass atoms.

set.seed(1)
state <- ripr_init(Q, plurality)
state <- sfw_step(state, times = 10L, until = gap_below(1e-10))
fit <- ripr_finish(
  state,
  correct_weights = TRUE,
  identify = TRUE,
  record_gap = TRUE
)

c(kl = fit@kl, gap = fit@gap_final, atoms = n_atoms(fit@W0))
#>         kl        gap      atoms 
#> 0.02822068 0.00000000 2.00000000

Plurality null, with the support of W1 (alternative) and W0 (null, fitted).

The fit attains a KL divergence (from Q) of approximately 0.0282. The gap quantity is a Frank–Wolfe duality gap: at its true value it bounds the expectation of Q / P* over the null by 1 + gap. The gap recorded during a fit (state@trace$gap_after, or fit@gap_fit) is only an estimate, so it is not itself an upper bound and cannot be used directly to get an e-variable. fit@gap_final, requested here via record_gap = TRUE, is a fresh sweep over the returned mixture, and is what we compare against the true certified bound below.

From a fit to an e-variable

With a fitted projection in hand, we can construct a random_variable that computes the mixture likelihood ratio X = Q / P*. This is our candidate e-variable:

X <- likelihood(Q, label = "Q") /
  likelihood(fit@P_star, label = "P*")
X
#> <random_variable> Q / P* 
#>   on count_space, dimension 3

random_variables are callable: they map elements of the sample space to their realisation.

outcomes <- rbind(
  c(10L, 10L, 0L),
  c(8L, 7L, 5L)
)
X(outcomes)
#> [1] 0.8372139 1.0794622

certify() computes an upper bound on the null expectation of X. It returns the upper bound alongside the largest value its internal search actually attained; the interval formed contains the true supremum, which is generally unknown.

cert <- certify(X, plurality, tol = 1e-9)
c(
  upper = cert@sup_ub,
  attained = cert@sup_lb,
  width = cert@sup_ub - cert@sup_lb
)
#>        upper     attained        width 
#> 1.000000e+00 1.000000e+00 9.455683e-10

The bound lands just above one: it is exactly 1 + gap, where gap is the true Frank–Wolfe duality gap that the fit stopped on.

c(bound_minus_one = cert@sup_ub - 1, fitted_gap = fit@gap_final)
#> bound_minus_one      fitted_gap 
#>    2.816769e-10    0.000000e+00

The rescaled random variable E = X / cert@sup_ub is then a genuine e-variable for plurality. e_variable() builds it from the certificate, dividing by the bound only when the bound is above one (a bound of one or less says X is an e-variable already):

E <- e_variable(cert)
print(E)
#> <random_variable> Q / P* / 1 
#>   on count_space, dimension 3
print(E(outcomes))
#> [1] 0.8372139 1.0794622

The correction doesn’t need a good fit

The whole point of splitting fitting from certification is that validity does not depend on the quality of the fit. To demonstrate, replace the projection with something that is not the RIPr at all (a uniform point mass) and certify the resulting ratio:

bad <- family(rep(1 / K, K))
X_bad <- likelihood(Q, label = "Q") /
  likelihood(bad, label = "P0")
cert_bad <- certify(X_bad, plurality, tol = 1e-9)

c(good = cert@sup_ub, bad = cert_bad@sup_ub)
#>     good      bad 
#>  1.00000 10.54509

Both rescale to valid e-variables. What differentiates them is e-power: the log-growth, or expected log e-value under the alternative:

growth <- function(E, Q) {
  outcomes <- enumerate_space(E@sample_space)
  e_values <- E(outcomes)
  Q_expectation <- log_density(Q, outcomes) |>
    exp() %*%
    log(e_values)
  c(Q_expectation) # return as scalar
}

E_gr <- growth(E, Q)

c(good = E_gr, bad = growth(e_variable(cert_bad), Q))
#>        good         bad 
#>  0.02822068 -1.99396747

A negative growth rate means evidence is lost on average even when the alternative is true. This means our e-variable is useless, but not invalid. Validity is free from the certification and correction procedure; fitting a good approximation is what buys e-power.

When P* is the true RIPr of Q onto the null, the attained KL divergence KL(Q || P*) is precisely the growth rate of the resulting e-variable. For any fitted P the KL divergence sits above that optimum, so fit@kl upper bounds the optimal growth rate. The true optimum therefore lies somewhere between the growth this fit achieved and fit@kl:

c(
  lower = E_gr,
  upper = fit@kl,
  width = fit@kl - E_gr
)
#>        lower        upper        width 
#> 2.822068e-02 2.822068e-02 2.816769e-10

A second worked example: composing step verbs

Because fw_step(), sfw_step(), lb_step(), em_step(), weight_step() and ct_step() are independent verbs rather than a fixed pipeline, they compose freely. Here we demonstrate a null where candidate 1, 2 and 3 all fail to reach a majority (H0: theta_1, theta_2, theta_3 <= 1/2, the medial triangle), fit with a loop that interleaves Frank–Wolfe and EM steps, and nominate an alternative with one dirac per section, equally weighted.

medial <- null_model(
  family,
  simplex_region(
    vertices = rbind(c(0.5, 0.5, 0), c(0.5, 0, 0.5), c(0, 0.5, 0.5))
  )
)

W_maj <- finite_dist(
  atoms = rbind(
    c(0.6, 0.2, 0.2),
    c(0.2, 0.6, 0.2),
    c(0.2, 0.2, 0.6)
  ),
  weights = c(1 / 3, 1 / 3, 1 / 3)
)
Q_maj <- family(W_maj)

Medial-triangle null, with the support of W_maj.

set.seed(2)
fit_medial <- ripr_init(Q_maj, medial) |>
  sfw_step(times = 40L, until = gap_below(1e-10)) |>
  ripr_finish(record_gap = TRUE)
X_medial <- likelihood(Q_maj, label = "Q") /
  likelihood(fit_medial@P_star, label = "P*")

c(kl = fit_medial@kl, upper = certify(X_medial, medial, tol = 1e-9)@sup_ub)
#>        kl     upper 
#> 0.3845465 1.0000000

Medial-triangle null, with the support of W_maj and W0 (fitted).

Searching vs. proving

sup_lb() hunts for the largest null expectation by multi-start gradient ascent. It is much cheaper than proving a bound, and it works on families and geometries that have no implemented bounding procedure, but it is not a gap certificate. It reports a lower bound on the supremum: something larger may sit somewhere gradient ascent could not reach.

c(searched = sup_lb(X, plurality)@sup_lb, certified = cert@sup_ub)
#>  searched certified 
#>         1         1

Where no bounding method has been implemented for a given combination of family and region geometry, certify() refuses outright:

gaussian <- gaussian_family(d = 2L)
null <- null_model(
  gaussian,
  halfspace_region(normal = c(1, -1), offset = 0)
)
Y <- likelihood(gaussian(c(0, 0)))
certify(Y, null)
#> Error:
#> ! Cannot certify:
#> No bounding method is implemented for gaussian_family expectations over halfspace_region.
#> Certifying this requires deriving and implementing a bound on gaussian_family expectations over halfspace_region. Nothing here says one does not exist. In the meantime, `sup_lb()` still searches, and reports a lower bound.

sup_lb() still works here as it needs no bounding method, only a search. The Gaussian sample space cannot be enumerated, so the expectation at each parameter is taken by quadrature instead: the same engine argument as ripr_init() takes, integrating under P_theta rather than under the alternative.

found <- sup_lb(
  Y,
  null,
  engine = gh_engine(15),
  n_seeds = 50L,
  n_restarts = 5L
)
found
#> <ripr_search>
#>   X = P[theta = (0, 0)]
#>   under null_model: gaussian_family over 1 part
#>   sup E[X] >= 0.07957749  (searched, not certified)
#>   attained at theta = (-2.83e-11, -4.5e-12), in part 1

Where to go next

The function reference lists every exported function and S7 class, with a runnable example on each. A few starting points if you’re extending rather than just using the package:

• ?parametric_family describes the interface a new observation model must implement: compile_loglik(), score() and kernel_draw(). ?sample_space and ?convex_region describe the two sets it maps between, and ?distribution the laws fam(theta) returns.
• vignette("regions") tours the region interface and the set algebra (|, &, -, == and disjoin()) used to build a null out of convex pieces. ?convex_region and its generics (chart(), contains(), project()) describe the interface a new convex geometry must implement; ?simplex_region and ?halfspace_region are the two built-in examples covered there.
• ?quadrature, ?exact_engine, ?mc_engine and ?gh_engine describe how expectations under the alternative are approximated, and how to add a new rule.
• ?oracles explains the Frank–Wolfe and Li–Barron oracles that sfw_step(), fw_step() and lb_step() are built on, and when to prefer one over the other.

References

Grünwald, Peter, Rianne de Heide, and Wouter Koolen. 2024. “Safe Testing.” Journal of the Royal Statistical Society Series B: Statistical Methodology 86 (5): 1091–128.

Li, Jonathan, and Andrew Barron. 1999. “Mixture Density Estimation.” Advances in Neural Information Processing Systems 12.