Computed, not simulated

Every number BEAM returns is computed from the model over every possible sequence, not estimated from simulated samples, so the same input always gives the same answer. The answers come in three kinds, and BEAM always tells you which kind you have.

Computed directly

Probabilities of sequences, the normaliser, length probabilities, and every sum over all sequences that a sweep through the model’s graph can compute directly; draws from the model follow its probabilities precisely. These carry only ordinary floating-point rounding, about one part in $10^{15}$ : the total probability of all sequences is one to that precision (see How BEAM works).

Computed with a guaranteed error bound

Quantities such as the expected number of distinct sequences at a depth, overlaps between samples and sequences private to one person are integrals that BEAM evaluates numerically. Each result comes as a bounded value: the number, together with a rigorous upper bound on its relative error. The bound is not an estimate of the error; the true value is guaranteed to lie within it.

from BEAM import BeamModel
from BEAM.analytics.exact_occupancy import ExactOccupancy

occupancy = ExactOccupancy(model)
r = occupancy.richness(50_000)
# the value, its guaranteed bound, and how it was computed
print(r)
print(f"{r.value:,.1f} distinct sequences, relative error bound {r.bound:.2g}")

try:
    # refused: the bound must not be dropped by accident
    float(r)
except TypeError as err:
    print("float():", err)

BoundedValue(49570.469422322, relative error <= 1.6e-12, contour q=0.78 h=0.03734)
49,570.5 distinct sequences, relative error bound 1.6e-12
float(): a BoundedValue is not a float; read .value (and .bound)

Full program, including the set-up hidden here

A bounded value cannot be turned into a plain number by accident: read .value once you have looked at .bound. Typical bounds lie between about $10^{-13}$ and $10^{-11}$ ; a value that could not be brought within the requested tolerance is returned with the bound it has, never silently.

Bounds instead of an approximation

A few quantities have no closed form for any model. For these BEAM returns what can be computed, guaranteed lower and upper bounds, rather than an approximation that could be mistaken for the value itself. The diversity of a group of repertoires is one: at Hill order 2 it is computed directly, but at order 1 (Shannon diversity) the logarithm of a sum of models does not break into per-position pieces, so BEAM brackets it instead:

from BEAM import BeamModel, Mixture

group = Mixture([donor_a, donor_b])            # pick a donor, then a cell
print(f"Hill number q=2: {group.hill_number(2.0):.6g}")

try:
    group.hill_number(1.0)                     # no closed form at q = 1
except Exception as err:
    print(f"q=1 declined: {type(err).__name__}")

bounds = group.entropy_bounds()                # guaranteed bounds instead
print(bounds)

Hill number q=2: 2.15324e+06
q=1 declined: UnsupportedOrder
EntropyBounds(lower=20.884631461318804, upper=21.577778641878748, weight_entropy=0.6931471805599453, exact=True)

Full program, including the set-up hidden here

The bracket is never wider than the entropy of the mixing weights, $\ln 2$ for two equal members.

What the computation covers

Computing what the model says, over every possible sequence, is one question. Whether the model describes your repertoire well is another: it depends on the data the model was built from and on the sampling assumptions (independent draws, no clonal expansion). The guides state, for each calculation, what it assumes.