cointf.naive: reference implementations¶
Reference (slow) exact likelihood for the Housworth-Stahl model.
A direct, one-meiosis-at-a-time implementation used to validate
cointf.fastlik. It uses the same generative model as
cointf.simintf:
Interference pathway: a stationary renewal process with gap density \(q f^*(qx; \nu)\), where \(q = 1 - p\) and
\[f^*(x; \nu) = \sum_{k \ge 1} 2^{-k}\, \mathrm{Gamma}(x;\, k\nu,\ \mathrm{rate}=2\nu).\]Escape pathway: a homogeneous Poisson process with rate \(p\).
Observed crossovers are the superposition of the two, with pathway labels unobserved. The likelihood sums over all labelings with a forward recursion on the position of the most recent interference crossover (\(O(n^2)\) per meiosis). \(f^*\), its survival function and its integrated tail are evaluated in log space from a 200-term gamma series.
loglik_quartet_naive() is the matching reference for phase-unknown
(quartet) data, the pooled crossovers of two meioses. It sums directly over
all \(3^n\) labelings (interference in meiosis A, interference in
meiosis B, or escape) and is used to validate cointf.quartet.
- cointf.naive.log_fstar(x, nu, kmax=200)¶
Log density of the unit-rate inter-crossover distance.
- Parameters:
- xarray_like
Distances in Morgans,
> 0.- nufloat
Interference parameter,
> 0.- kmaxint, default 200
Number of gamma-series terms.
- Returns:
- numpy.ndarray
\(\log f^*(x; \nu)\), shape
(x.size,).
- cointf.naive.log_sstar(x, nu, kmax=200)¶
Log survival function of the inter-crossover distance.
- Parameters:
- xarray_like
Distances in Morgans,
>= 0.- nufloat
Interference parameter,
> 0.- kmaxint, default 200
Number of gamma-series terms.
- Returns:
- numpy.ndarray
\(\log(1 - F^*(x; \nu))\), shape
(x.size,).
- cointf.naive.log_tail_integral(x, nu, kmax=200)¶
Log integrated survival function of the inter-crossover distance.
Computes \(\log G(x)\) with \(G(x) = \int_x^\infty (1 - F^*(u; \nu))\, du\). For a gamma component with shape \(a\) and rate \(b\), \(\int_x^\infty S(u)\, du = (a/b) Q(a + 1, bx) - x Q(a, bx)\), with \(Q\) the regularised upper incomplete gamma function; \(G(0) = 1\).
- Parameters:
- xfloat
Length in Morgans,
>= 0.- nufloat
Interference parameter,
> 0.- kmaxint, default 200
Number of gamma-series terms.
- Returns:
- float
\(\log G(x)\), or
-infif the series is not positive.
- cointf.naive.loglik_meiosis_naive(xs, L, nu, p, kmax=200)¶
Exact log-likelihood of one meiosis’ crossover positions.
- Parameters:
- xsarray_like
Crossover positions in Morgans, in
[0, L]and distinct.- Lfloat
Chromosome length in Morgans,
> 0.- nufloat
Interference parameter,
> 0.- pfloat
Fraction of crossovers escaping interference, in
[0, 1].- kmaxint, default 200
Number of gamma-series terms.
- Returns:
- float
Log of the joint density of the observed crossover positions.
- Raises:
- ValueError
If L, nu or p is out of range.
See also
cointf.fastlik.loglik_fastFast, vectorised equivalent.
- cointf.naive.loglik_naive(data, L, nu, p, kmax=200)¶
Exact log-likelihood summed over independent meioses.
- Parameters:
- datasequence of array_like
Crossover positions in Morgans, one sequence per meiosis.
- Lfloat
Chromosome length in Morgans, shared by all meioses.
- nufloat
Interference parameter,
> 0.- pfloat
Fraction of crossovers escaping interference, in
[0, 1].- kmaxint, default 200
Number of gamma-series terms.
- Returns:
- float
Sum of
loglik_meiosis_naive()over the meioses.
- cointf.naive.loglik_quartet_naive(xs, L, nu, p, kmax=200)¶
Brute-force log-likelihood of one phase-unknown (quartet) dataset.
Sums directly over all \(3^n\) labelings of the observed crossovers as interference crossovers of meiosis A, of meiosis B, or escape crossovers (a rate-\(2p\) Poisson process). Only practical for small \(n\); used to validate
cointf.quartet.loglik_quartet().- Parameters:
- xsarray_like
Pooled crossover positions of two meioses, in
[0, L]and distinct.- Lfloat
Chromosome length in Morgans,
> 0.- nufloat
Interference parameter,
> 0.- pfloat
Escape fraction in each meiosis, in
[0, 1].- kmaxint, default 200
Number of gamma-series terms.
- Returns:
- float
Log of the joint density of the pooled crossover positions.
- Raises:
- ValueError
If L, nu or p is out of range.