cointf.fastlik: single-meiosis likelihood¶

Fast exact likelihood for the Housworth-Stahl model.

Same model and recursion as cointf.naive, organised for speed:

\(f^*\), \(S^* = 1 - F^*\) and the integrated tail \(G\) are gamma series. For each distance only the terms that matter are evaluated: an index window is found by bisection with explicit error bounds, and the terms above it are summed in closed form. Typically about 10 of the terms need an incomplete gamma or log-gamma evaluation.
The recursion needs \(\log f^*\) only at gaps between crossovers and \(\log S^*\) only at the distances to the chromosome ends, so each is evaluated only where it is used.
Meioses are grouped by crossover count and the forward recursion runs once per group on a dense array, so cost is a few numpy calls per crossover index and memory scales with the meioses present.
All kernel evaluations needed for one likelihood call are made in a single vectorised call, so the kernel can be swapped without touching the recursion.

The interference parameter \(\nu\) and escape fraction \(p\) may be scalars (a pooled fit) or arrays with one value per meiosis (e.g. when they depend on covariates).

Notes¶

Interference crossovers form a stationary renewal process with gap density \(q f^*(qx; \nu)\), \(q = 1 - p\), where

\[f^*(x; \nu) = \sum_{k \ge 1} 2^{-k}\, \mathrm{Gamma}(x;\, k\nu,\ \mathrm{rate}=2\nu).\]

Escape crossovers form an independent Poisson process of rate \(p\) (per Morgan). Observed crossovers are the superposition, with pathway labels unobserved; the likelihood sums over all labelings with an \(O(n^2)\) forward recursion per meiosis with \(n\) crossovers.

Truncation. With \(z = 2\nu x\) and shape \(a = k\nu\), the series terms are bounded with the Chernoff exponent \(e(a, z) = z - a + a \log(a / z)\): \(Q(a, z) \le e^{-e}\) for \(a < z\) and \(1 - Q(a, z) \le e^{-e}\) for \(a > z\), where \(Q\) is the regularised upper incomplete gamma function. Terms are dropped, or replaced by their limit \(Q = 1\), only where the total error is below \(e^{-40}\) relative to a lower bound on the result. The log terms of \(f^*\) are concave in \(k\), so its window is the set of terms within a fixed margin of the largest one. Series indices are capped at max_series_terms.

class cointf.fastlik.Meioses(data, L)¶

Crossover data prepared for repeated likelihood evaluation.

Validates and sorts the positions, groups meioses by crossover count and precomputes every distance the likelihood needs. Build it once and pass it to loglik_fast() for each parameter value.

Parameters:

datasequence of array_like: Crossover positions in Morgans, one 1-D sequence per meiosis (a meiosis here is one chromosome in one gamete). Positions must be finite, lie in [0, L] and be distinct within a meiosis; order does not matter. For cointf.quartet.loglik_quartet() each entry is instead one quartet, the pooled crossovers of two meioses of the same chromosome, and “meiosis” in the attributes below means one such entry.
Lfloat or array_like: Chromosome length in Morgans, > 0: a scalar shared by all meioses or one value per meiosis.

Attributes:

Mint: Number of meioses.
Lnumpy.ndarray: Chromosome length per meiosis, shape (M,).
countsnumpy.ndarray: Number of crossovers per meiosis, shape (M,).
groupslist of tuple: One (n, idx, iu, ju, end_bounds, gap_bounds) per crossover count n >= 1: meiosis indices, upper-triangle pair indices, and slice bounds into x_end and x_gap.
x_endnumpy.ndarray: Distances from each crossover to the chromosome start and end, where \(S^*\) is needed, before scaling by \(q\).
row_endnumpy.ndarray: Meiosis index of each entry of x_end.
x_gapnumpy.ndarray: Distances between pairs of crossovers, where \(f^*\) is needed, before scaling by \(q\).
row_gapnumpy.ndarray: Meiosis index of each entry of x_gap.
no_crossovernumpy.ndarray: Indices of meioses without crossovers.

Raises:

ValueError: If L has the wrong shape or is not finite and positive, or a meiosis has non-finite, out-of-range or duplicated positions.

Examples

>>> from cointf.fastlik import Meioses, loglik_fast
>>> mz = Meioses([[0.2, 0.9], [], [0.5]], L=1.2)
>>> mz.M, mz.counts.tolist()
(3, [2, 0, 1])
>>> ll = loglik_fast(mz, None, nu=4.0, p=0.05)

cointf.fastlik.log_fstar_fast(x, nu)¶

Log density of the unit-rate inter-crossover distance.

Parameters:

xarray_like: Distances in Morgans, >= 0. Flattened to 1-D.
nufloat or array_like: Interference parameter, > 0. A scalar, or an array with the same number of elements as x (one value per distance).

Returns:

numpy.ndarray: \(\log f^*(x; \nu)\), shape (x.size,). At x = 0 this is -inf for nu > 1, 0 for nu = 1 and +inf for nu < 1.

See also

cointf.naive.log_fstar: Slower 200-term reference implementation.

Examples

>>> from cointf.fastlik import log_fstar_fast
>>> log_fstar_fast([0.5, 1.0], nu=1.0)  # nu = 1 is Exp(1)
array([-0.5, -1. ])

cointf.fastlik.log_sstar_fast(x, nu)¶

Log survival function of the inter-crossover distance.

Evaluates \(\log S^*(x; \nu) = \log(1 - F^*(x; \nu))\).

Parameters:

xarray_like: Distances in Morgans, >= 0. Flattened to 1-D.
nufloat or array_like: Interference parameter, > 0. A scalar, or an array with the same number of elements as x.

Returns:

numpy.ndarray: \(\log S^*(x; \nu)\), shape (x.size,); 0 at x = 0.

See also

cointf.naive.log_sstar: Slower 200-term reference implementation.

cointf.fastlik.log_tail_fast(x, nu)¶

Log integrated survival function of the inter-crossover distance.

Evaluates \(\log G(x; \nu)\) with \(G(x) = \int_x^\infty S^*(u; \nu)\, du\), using the closed form \(\int_x^\infty S_k(u)\, du = (a/b) Q(a + 1, bx) - x Q(a, bx)\) for each gamma component (shape \(a\), rate \(b\), \(Q\) the regularised upper incomplete gamma function). \(G(0) = 1\).

Parameters:

xarray_like: Lengths in Morgans, >= 0. Flattened to 1-D.
nufloat or array_like: Interference parameter, > 0. A scalar, or an array with the same number of elements as x.

Returns:

numpy.ndarray: \(\log G(x; \nu)\), shape (x.size,). -inf where the series underflows to a non-positive value.

See also

cointf.naive.log_tail_integral: Scalar 200-term reference implementation.

cointf.fastlik.loglik_fast(data, L, nu, p, *, kernel=None, per_meiosis=False)¶

Exact log-likelihood of crossover positions under the Housworth-Stahl model.

Parameters:

datasequence of array_like or Meioses: Crossover positions in Morgans, one sequence per meiosis, or a prepared Meioses object (faster for repeated evaluation).
Lfloat or array_like or None: Chromosome length(s) in Morgans: a scalar or one value per meiosis. Must be None when data is a Meioses, which carries its own lengths.
nufloat or array_like: Interference parameter, > 0: a scalar shared by all meioses or one value per meiosis. nu = 1 is no interference.
pfloat or array_like: Fraction of crossovers escaping interference, in [0, 1]: a scalar or one value per meiosis.
kernelcallable, optional: Function kernel(u_gap, nu_gap, u_end, nu_end) -> (log_f, log_S) returning \(\log f^*\) at the gaps u_gap and \(\log S^*\) at the end distances u_end (see series_kernel(), the default, and cointf.chebyshev.ChebyshevKernel, a faster approximation). The nu arguments are scalars when nu is shared by all meioses, otherwise arrays matching their distances.
per_meiosisbool, default False: If True, return one log-likelihood per meiosis instead of the sum.

Returns:

float or numpy.ndarray: The total log-likelihood, or an array of shape (M,) if per_meiosis is True.

Raises:

ValueError: If parameters are out of range or have the wrong shape, if L is given together with a Meioses, or if the data are invalid (see Meioses).

See also

Meioses: Prepare data once for repeated evaluation.
cointf.naive.loglik_meiosis_naive: Slow single-meiosis reference.

Notes

Passing an array whose entries are all equal is treated as the corresponding scalar, so a pooled fit costs the same either way.

Examples

>>> from cointf.fastlik import loglik_fast
>>> data = [[0.31, 1.12], [0.70], []]
>>> total = loglik_fast(data, 1.5, nu=4.0, p=0.05)
>>> per = loglik_fast(data, 1.5, nu=[4.0, 6.0, 2.0], p=0.05, per_meiosis=True)
>>> per.shape
(3,)

cointf.fastlik.series_kernel(u_gap, nu_gap, u_end, nu_end)¶

Default likelihood kernel: exact \(\log f^*\) and \(\log S^*\).

The recursion needs \(f^*\) only at gaps between interference crossovers and \(S^*\) only at distances to the chromosome ends, so a kernel receives the two sets separately. Any callable with this signature can be passed as kernel to loglik_fast() or cointf.quartet.loglik_quartet(). cointf.chebyshev.ChebyshevKernel is a faster, approximate alternative.

Parameters:

u_gapnumpy.ndarray: 1-D array of gaps between crossovers on the interference scale (Morgans times \(q = 1 - p\)), >= 0.
nu_gapfloat or numpy.ndarray: Interference parameter: a scalar, or an array shaped like u_gap.
u_endnumpy.ndarray: 1-D array of distances from a crossover to a chromosome end, on the interference scale, >= 0.
nu_endfloat or numpy.ndarray: Interference parameter: a scalar, or an array shaped like u_end.

Returns:

log_fnumpy.ndarray: \(\log f^*(u; \nu)\) at u_gap.
log_Snumpy.ndarray: \(\log S^*(u; \nu)\) at u_end.