Introduction
A null hypothesis, or null_model in ripr
comprises a parametric_family and a region: a
subset of the family’s parameter_space.
convex_regions are individual convex units that carry
the geometry a fitting procedure needs: a chart() provides
a coordinate system for optimisation, contains() encodes
membership for the set, and project() moves outside points
onto the region. The union_region is a finite union of
convex regions and is generally not convex.
s <- simplex_region(vertices = diag(3))
S7::S7_inherits(s, space)
#> [1] TRUE
S7::S7_inherits(s, region)
#> [1] TRUE
S7::S7_inherits(s, convex_region)
#> [1] TRUE
# a space with no geometry is not a region.
S7::S7_inherits(count_space(n_trials = 4L, k = 3L), space)
#> [1] TRUE
S7::S7_inherits(count_space(n_trials = 4L, k = 3L), region)
#> [1] FALSE
u <- union_region(s, halfspace_region(normal = c(1, -1, 0)))
S7::S7_inherits(u, region)
#> [1] TRUE
S7::S7_inherits(u, convex_region)
#> [1] FALSEThere are two methods that all unions support: parts()
and cells(). parts() returns all the convex
pieces that the region was declared as, whereas
cells() returns a triangulation. The two agree for a
union_region, but not for a polytope_region
that isn’t already a simplex: the certification routine needs simplices
while fitting can operate directly on polytopes, so a square region has
one part and two triangular cells, each filling different roles in the
package.
Convex geometries
Every convex region in the package is, underneath, a
polyhedron_region(): the Minkowski–Weyl generator form
conv(V) + cone(R) + span(L) (or its dual facet description,
if you build one from h_region() instead).
simplex_region(), polytope_region(),
point_region(), real_region() and
halfspace_region() special cases of convex polyhedra, we
provide them as convenient constructors for more familiar objects, with
validation. simplex_region() and
halfspace_region() are the two you’ve already seen in
vignette("ripr"). Here are the rest:
# The convex hull of a vertex set, one vertex per row. A simplex_region is
# the special case when the vertices are affinely independent.
polytope_region(vertices = rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1)))
#> <polytope_region>
#> 4 vertices in R^2
#> facets: 4 inequalities
#> vertices:
#> [,1] [,2]
#> [1,] 0 0
#> [2,] 1 0
#> [3,] 1 1
#> [4,] 0 1
# A single parameter point, e.g. to specify a simple null theta = (0.5 0.3 0.2)'
point_region(theta = c(0.5, 0.3, 0.2))
#> <point_region>
#> the point (0.5, 0.3, 0.2)
#> facets: 3 equalities
# All of R^2
real_region(2L)
#> <real_region>
#> all of R^2
#> facets: noneSet algebra
regions are composed with the operators |
(union), & (intersection), - (difference)
and == (equality), along with disjoin(), all
working in exact rational arithmetic so overlaps and boundaries are
decided precisely rather than approximately. Base R’s
union(), intersect(), setdiff()
and setequal() are left untouched, and do not accept
regions. The plurality null from vignette("ripr") was
constructed from two overlapping triangles.
K <- 3L
plurality_part <- function(j) {
vertices <- diag(K)
vertices[1L, ] <- replace(numeric(K), c(1L, j), 0.5)
simplex_region(vertices = vertices)
}
plurality <- plurality_part(2) | plurality_part(3)
plurality
#> <union_region>
#> 2 parts, dimension 3
#> simplex_region: 3 vertices in R^3
#> simplex_region: 3 vertices in R^3Taking the intersection (via &) instead yields the
region where the two parts overlap: the wedge where
theta_1 <= theta_2, theta_3.
plurality_part(2) & plurality_part(3)
#> <simplex_region>
#> 3 vertices in R^3
#> facets: 3 inequalities, 1 equality
#> vertices:
#> [,1] [,2] [,3]
#> [1,] 0.0000000 0.0000000 1.0000000
#> [2,] 0.3333333 0.3333333 0.3333333
#> [3,] 0.0000000 1.0000000 0.0000000x - y is the closure of x minus
y. The closed complement of the plurality null within the
full simplex is the region where
theta_1 >= theta_2, theta_3. It is not quite the
alternative because it includes the boundary where
theta_1 = theta_2 or theta_1 = theta_3.
leader <- simplex_region(vertices = diag(K)) - plurality
leader
#> <polytope_region>
#> 4 vertices in R^3
#> facets: 4 inequalities, 1 equality
#> vertices:
#> [,1] [,2] [,3]
#> [1,] 0.5000000 0.0000000 0.5000000
#> [2,] 0.3333333 0.3333333 0.3333333
#> [3,] 0.5000000 0.5000000 0.0000000
#> [4,] 1.0000000 0.0000000 0.0000000The plurality null’s two parts overlap, which is fine for fitting and
certification (a supremum over a union is a supremum over any cover of
it), but sometimes you may wish for a disjoint cover instead. For
instance, to sum a measure over the null without double-counting the
overlap. disjoin() peels a union apart sequentially: it
leaves the first part as is, then subtracts the first part from the
second, and so on.
peeled <- disjoin(plurality)
length(parts(plurality))
#> [1] 2
length(parts(peeled))
#> [1] 2
parts(plurality)[[2]]
#> <simplex_region>
#> 3 vertices in R^3
#> facets: 3 inequalities, 1 equality
#> vertices:
#> [,1] [,2] [,3]
#> [1,] 0.5 0 0.5
#> [2,] 0.0 1 0.0
#> [3,] 0.0 0 1.0
parts(peeled)[[2]]
#> <simplex_region>
#> 3 vertices in R^3
#> facets: 3 inequalities, 1 equality
#> vertices:
#> [,1] [,2] [,3]
#> [1,] 0.0000000 0.0000000 1.0000000
#> [2,] 0.3333333 0.3333333 0.3333333
#> [3,] 0.5000000 0.0000000 0.5000000== decides whether two regions describe the same set,
returning a single TRUE or FALSE. This is
computed exactly, regardless of how each set is decomposed into parts,
so it compares sets rather than objects; identical() still
tells you whether two regions were built the same way.
peeled == plurality
#> [1] TRUE
identical(peeled, plurality)
#> [1] FALSE
square == union_region(cells(square))
#> [1] TRUE
square == cells(square)[[1]]
#> [1] FALSEWhen nothing is left, the algebra says so rather than erroring out.
The intersection of disjoint regions, or a difference that removes
everything, is an empty_region(): it behaves like a list of
length zero, and | drops it without a trace.
nothing <- point_region(theta = c(1, 0, 0)) & point_region(theta = c(0, 1, 0))
is_empty(nothing)
#> [1] TRUE
identical(nothing | s, s)
#> [1] TRUEA Sierpinski-like null
vignette("ripr") builds the medial null for a
three-category multinomial: the null region is where no candidate holds
a majority, H0: theta_1, theta_2, theta_3 <= 1/2.
Geometrically it is the medial triangle of the 2-simplex, the one you
get by joining the midpoints of each edge. The alternative is the union
of the three corner triangles.
The corners themselves are easily obtained by taking the complement
of the null within the full simplex with -:
medial_of <- function(triangle) {
v <- triangle@vertices
simplex_region(
vertices = rbind(
(v[1, ] + v[2, ]) / 2,
(v[1, ] + v[3, ]) / 2,
(v[2, ] + v[3, ]) / 2
)
)
}
full <- simplex_region(vertices = diag(3))
medial <- medial_of(full)
corners <- parts(full - medial)
length(corners)
#> [1] 3Repeating the same split inside each corner (take its own medial
triangle, add that to the null) a few times yields a Sierpinski
approximation built entirely from | and -:
sierpinski_null <- function(triangle, depth) {
medial <- medial_of(triangle)
if (depth <= 1L) {
return(medial)
}
corners <- parts(triangle - medial)
medial | union_region(lapply(corners, sierpinski_null, depth = depth - 1L))
}
depth <- 3L
geometry <- sierpinski_null(full, depth)
length(parts(geometry))
#> [1] 13The finest, still-uncovered corners are the natural place to put an
alternative – one dirac per leftover triangle, exactly as
vignette("ripr") puts one per corner of the plain medial
null, just at finer resolution:
leaf_corners <- function(triangle, depth) {
if (depth <= 0L) {
return(list(triangle))
}
corners <- parts(triangle - medial_of(triangle))
unlist(lapply(corners, leaf_corners, depth = depth - 1L), recursive = FALSE)
}
leaves <- leaf_corners(full, depth)
length(leaves)
#> [1] 27
centroids <- do.call(rbind, lapply(leaves, \(cell) colMeans(cell@vertices)))
W_alt <- finite_dist(
atoms = centroids,
weights = rep(1 / length(leaves), length(leaves))
)
A depth-3 Sierpinski-like null, with the support of the alternative – one dirac per leftover corner.
Every part of geometry is still just a
simplex_region, so nothing about fitting or certification
needs to know it came from a recursive construction rather than being
typed out by hand as two parts:
set.seed(3)
fit <- ripr_init(Q_alt, sierpinski) |>
sfw_step(times = 40L, until = gap_below(1e-10)) |>
ripr_finish(correct_weights = TRUE, identify = TRUE, record_gap = TRUE)
c(kl = fit@kl, atoms = n_atoms(fit@W0))
#> kl atoms
#> 0.001768473 42.000000000
The same null, with the support of the alternative and the fitted null mixture W0.
X <- likelihood(Q_alt, label = "Q") /
likelihood(fit@P_star, label = "P*")
cert <- certify(X, sierpinski, tol = 1e-9)
c(upper = cert@sup_ub, attained = cert@sup_lb)
#> upper attained
#> 1 1certify() proves the same kind of bound it did for the
two-part plurality null in vignette("ripr"), over a
union_region of 13 parts built by a few lines of recursive
set algebra rather than declared by hand.
E <- e_variable(cert) is a genuine e-variable for
sierpinski, exactly as before. The gap is quite large,
however, which indicates that it is probably quite a poor fit.
Where next
vignette("ripr") covers the other half: fitting a
projection, why validity survives a deliberately bad one, and what
certify() refuses to do without a bounding method. If
you’re implementing a new convex geometry rather than composing the
built-in ones, ?convex_region and its generics
(chart(), contains(), project())
describe the interface it must satisfy.