POSQ

The content on this page was written by AI under human supervision.

POSQ computes the Bayesian evidence of a small continuous-time Markov chain model on a tree (the likelihood integrated over the prior) and returns an interval whose lower and upper endpoints are both proven bounds. It reduces the integral to an exact Gauss quadrature sum in which every term is positive and evaluates that sum in ball arithmeticevery number is carried as a midpoint plus a radius that provably contains the true value; FLINT/Arb is the library used. The inputs are site-pattern counts for a four-leaf tree and the prior's parameters; the outputs are a certified interval for the evidence of each tree shape and certified lower bounds on log Bayes factors between shapes.

What it does

Bayesian model comparison rests on the evidence $Z=\int(\text{likelihood})\times(\text{prior})$. For tree models it is usually estimated by Monte Carlo, which gives a number with a statistical error and no guaranteed bound on either side. Interval arithmetic applied directly gives uselessly wide bounds, because expanding the likelihood into polynomial coefficients produces enormous cancellation. POSQ avoids both by never expanding anything: the likelihood is evaluated as a product of values at quadrature nodes, and every quantity entering the sum is positive.

The model class is a two-state Markov chain on a rooted four-leaf tree with a strict clock, exponential priors (rate $\lambda$) on the clock increments $s_k$, and a uniform prior on the stationary frequency $p$. Substituting $u_k=e^{-\beta s_k}$, with $\beta$ the chain's rate constant, integrates the exponential priors out exactly and leaves

$$Z_s(p)=\frac{c^3}{2}\int_{[0,1]^3}u^{\,c-1}\,F(u;p)\,du,\qquad c=2\lambda p(1-p),$$

where $F$ is the product over observed site patterns of the pattern probability raised to its count: a polynomial in $u$ of known degree, positive on the cube. A Gauss–Jacobi rulea quadrature rule with n nodes that integrates polynomials up to degree 2n−1 exactly against the weight u^(c−1) on [0,1] with enough nodes integrates $F$ with zero error, so the three clock dimensions add nothing to the width of the answer. POSQ builds the nodes and weights by interval Newton iteration at 2048-bit precision and accepts a rule only when every weight is provably positive and the rule reproduces exact moments. With positive weights and positive values there is no cancellation: a full-size sweep of 48.7 million nodes at 192 bits encloses $Z_s(p)$ to relative width near $2^{-179}$ in about 0.2 CPU-hours.

The remaining dimension, $p$, is integrated over fixed panels with certified Gauss–Legendre nodes, a remainder bounded by an exactly integrable Bernstein-basis majorant, and exact bands at the endpoints. All of the final width and nearly all of the cost are here. Measured on a 318-site four-taxon data set, one tree shape costs 27 to 38 CPU-hours for an $\ln Z$ interval 0.06 to 0.1 nats wide, and a looser target barely lowers the cost. A certified comparison is then the winner's lower endpoint minus the rival's upper endpoint, rounded outward for display.

The code raises an error instead of returning a number when a weight ball straddles zero, when a per-node value or an assembled $Z$ fails to certify positive, or when a sweep's relative width is worse than $2^{-100}$. The three balanced shapes $((A,B),(C,D))$ use a separate min/difference change of variables whose exactness checks must pass before any sweep.

Limits. Rooted four-leaf trees only (12 caterpillar and 3 balanced shapes); six leaves and ascertainment-corrected likelihoods are not covered. A different model in the same family needs its degree bookkeeping and prior-to-weight map re-derived, its rules rebuilt, and the verification rerun. The repository holds the engine, an exact-rational 8-count surrogate problem for testing, and the verification driver. The full-size driver, data and result tables are not included; scripts that need the data stop with POSQ-PRODUCTION-DATA-ABSENT unless POSQ_STAGE0 points at a directory containing it. A separate module, closure_interface.py, estimates the evidence of a user-written log-likelihood on the unit cube by Monte Carlo; its numbers are floating-point estimates, not certified, and it needs an upstream build (POSQ_R3_ROOT) that is not in the repository.

Examples

Quick self-test. Check that the exact surrogate and the C kernel agree on this machine:

python3 posq.py --selftest

S1 recomputes the surrogate's evidence at $(p,\lambda)=(709/2048,\,10)$ as an exact fraction and compares its checksum and logarithm with the recorded values. S2 runs the kernel at the same point and requires the certified ball to contain that exact rational at relative width better than $2^{-100}$. A passing run prints

[PASS] S1 exact surrogate sha256 pin (e6fc044cdd0dabfb..)
[PASS] S1b ln Z = -25.97566273247789 vs pinned -25.97566273247789
[PASS] S2 live kernel ball contains exact rational (relwidth 2^-248.2, line 2^-100)
posq selftest: PASS

and exits 0. If no kernel binary runs here and none can be built, S2 is skipped with the message POSQ-KERNEL-UNAVAILABLE.

One certified value from Python. Evaluate the surrogate's evidence through the kernel at an exact rational point:

import sys; sys.path.insert(0, '<your-checkout>/tools/posq')
import posq

adapt = posq.PosqKernelAdaptation()
Z = adapt.point_eval(("709/2048", "10"))   # certified ball via the C kernel

Z is an Arb ball: the true value lies within Z.rad() of Z.mid(). Before sweeping, the object rebuilds the three Gauss–Jacobi rules at this point's exact $c$ and rechecks them; a rule that fails raises rather than returning a loose ball. adapt.deriv(p, dim) gives $\partial Z/\partial p$ or $\partial Z/\partial\lambda$ from an exact closed form with no quadrature, which is the independent path the verification compares against.

Full verification run. Put the engine through the ERAS adversarial verification suite:

python3 verify_posq_adaptation.py [--seed N] [--points N]

Phase 0 repeats the self-test, adds a deliberately corrupted count vector that the kernel must detect, and checks that a box enclosure contains the exact value at the box's center, corners and edge midpoints. Phase 1 samples parameter boxes on three shells and requires every enclosure to contain the kernel's point values. It then plants narrowed, shifted and dimension-dropped enclosures that must be caught, rebuilds the rules at 53 bits to confirm that too little precision yields no answer rather than a wrong one, and compares finite differences with the exact derivatives. Exit 0 is pass, 1 fail, 2 indeterminate (for example, no runnable kernel). One passing seed is necessary, not sufficient; rerun with several --seed values before trusting an adaptation of your own.

Routines

Command line

Python engine

Closure interface (closure_interface.py)

Requirements and source

Python 3 with python-flint; scipy for the initial node guesses in the rule builder; numpy and scipy for the closure-interface self-test. No kernel binary is included: on first use kernel_io compiles posq_kernel.c into a cache directory (POSQ_KERNEL_CACHE, default ~/.cache/posq), which needs a C compiler plus the FLINT, MPFR and GMP headers. You can also build it yourself with gcc -O2 -o posq_kernel posq_kernel.c -lflint -lmpfr -lgmp -lm and set POSQ_KERNEL. Self-tests: python3 posq.py --selftest and python3 verify_posq_adaptation.py (needs tools/eras/ beside it). The code is in tools/posq/ in BootLoops' bootloops-dev repository (GitHub organization BootLoops-ai), released under the MIT license, and is also importable as the quad.posq module of Baller, which links to the same files.

← back to the tools index