cointf.quartet: phase-unknown likelihood

Phase-unknown (quartet) likelihood for the Housworth-Stahl model.

In a nuclear family with two offspring and no grandparents, a parent’s crossovers are seen as switches in the siblings’ shared inheritance, but it is unknown in which of the two meioses each one happened. The observed positions are then the superposition of two independent meioses with the same \((\nu, p)\) and length \(L\).

Each observed crossover comes from the interference pathway of meiosis A, from that of meiosis B, or from the escape pathway of either meiosis. The two escape pathways superpose to a Poisson process of rate \(2p\), so escape crossovers need no meiosis label. The likelihood sums over all \(3^n\) labelings with a forward recursion whose state is the pair of latest interference crossovers in A and B, costing \(O(n^3)\) per quartet instead of \(O(3^n)\).

The interference terms are the same as for a single meiosis (\(f^*\) at gaps, \(S^*\) at end distances and \(G\) for an empty pathway), so this module reuses cointf.fastlik.Meioses and the kernel interface of cointf.fastlik.loglik_fast().

cointf.quartet.loglik_quartet(data, L, nu, p, *, kernel=None, per_meiosis=False)

Log-likelihood of phase-unknown crossovers pooled from two meioses.

Parameters:

datasequence of array_like or Meioses

Crossover positions in Morgans, one sequence per quartet (the pooled crossovers of two meioses of one chromosome), or a prepared cointf.fastlik.Meioses.

Lfloat or array_like or None

Chromosome length(s) in Morgans; None when data is a cointf.fastlik.Meioses.

nufloat or array_like

Interference parameter, > 0: a scalar or one value per quartet. Both meioses of a quartet share it.

pfloat or array_like

Escape fraction in each meiosis, in [0, 1]: a scalar or one value per quartet. Both meioses of a quartet share it.

kernelcallable, optional

Kernel kernel(u_gap, nu_gap, u_end, nu_end) -> (log_f, log_S); see cointf.fastlik.series_kernel() (the default) and cointf.chebyshev.ChebyshevKernel.

per_meiosisbool, default False

If True, return one log-likelihood per entry of data, i.e. per quartet (one pair of pooled meioses), instead of the sum. Named as in cointf.fastlik.loglik_fast() to keep the two interfaces the same; it never splits a quartet into its two meioses.

Returns:

float or numpy.ndarray

The total log-likelihood, or an array of shape (M,) with one value per quartet if per_meiosis is True, where M is the number of quartets.

Raises:

ValueError

As for cointf.fastlik.loglik_fast().

See also

cointf.naive.loglik_quartet_naive

Brute-force \(3^n\) reference.

cointf.simintf.simulate_quartets

Simulate quartet data.

Notes

With nu = 1 or p = 1 the data are a rate-2 Poisson process and the log-likelihood is \(n \log 2 - 2L\).

Examples

>>> from cointf.quartet import loglik_quartet
>>> data = [[0.2, 0.31, 0.9, 1.12], [0.7], []]
>>> ll = loglik_quartet(data, 1.5, nu=4.0, p=0.05)