Expected site frequency spectrum
Random joining on block sizes
DefinitionMerging two blocks of a block-size list#✓fully proved
Merge the blocks at positions i < j of the block-size list b: remove both and append a
block of size b[i] + b[j].
DefinitionUnordered pairs of blocks#✓fully proved
The unordered pairs of positions i < j among k blocks.
DefinitionRandom-joining law on block-size lists#✓fully proved
The random-joining law: randomJoining m b ps is the probability that m steps of
random joining started from the block-size list b visit the list of block-size lists ps
(starting with b itself). Each step merges a pair of blocks chosen uniformly among all
pairs; the process stops once fewer than two blocks remain.
noncomputable def randomJoining : ℕ → List ℕ → List (List ℕ) → ℝ
| 0, b, ps => if ps = [b] then 1 else 0
| m + 1, b, ps =>
if (blockPairs b.length).card = 0 then (if ps = [b] then 1 else 0)
else match ps with
| [] => 0
| b' :: rest =>
if b' = b then
(∑ p ∈ blockPairs b.length, randomJoining m (mergeBlocks b p.1 p.2) rest) /
(blockPairs b.length).card
else 0
def blockPairs (k : ℕ) : Finset (ℕ × ℕ) :=
(range k ×ˢ range k).filter fun p => p.1 < p.2
MathlibExt/Biology/PopulationGenetics/Coalescent/SiteFrequencySpectrum.lean:77
def mergeBlocks (b : List ℕ) (i j : ℕ) : List ℕ :=
((b.eraseIdx j).eraseIdx i) ++ [b.getD i 0 + b.getD j 0]
MathlibExt/Biology/PopulationGenetics/Coalescent/SiteFrequencySpectrum.lean:73
Uses 2 · Used by 4
DefinitionTransition operator of random joining#✓fully proved
The transition operator of random joining acting on a statistic g of block-size lists:
the average of g over the lists obtained from b by merging a uniformly chosen pair of
blocks (g [] once fewer than two blocks remain). By exists_finset_randomJoining, its j-th
iterate at b is the expectation of g at step j of random joining started from b.
noncomputable def joinStep (g : List ℕ → ℝ) (b : List ℕ) : ℝ :=
if (blockPairs b.length).card = 0 then g []
else (∑ p ∈ blockPairs b.length, g (mergeBlocks b p.1 p.2)) / (blockPairs b.length).card
def blockPairs (k : ℕ) : Finset (ℕ × ℕ) :=
(range k ×ˢ range k).filter fun p => p.1 < p.2
MathlibExt/Biology/PopulationGenetics/Coalescent/SiteFrequencySpectrum.lean:77
def mergeBlocks (b : List ℕ) (i j : ℕ) : List ℕ :=
((b.eraseIdx j).eraseIdx i) ++ [b.getD i 0 + b.getD j 0]
MathlibExt/Biology/PopulationGenetics/Coalescent/SiteFrequencySpectrum.lean:73
DefinitionTotal branch length subtending i sample copies#✓fully proved
The total length of the branches subtending exactly i of the n sample copies, for the
block-size history ps (entry n - k lists the block sizes while there are k lineages)
and holding times T: ∑_{k=2}^{n} #{blocks of size i among k} · T k.
Finite support and marginals of the random-joining law
TheoremRandom joining has finite support and iterates the join step#✓fully proved
The random-joining law has finite support, and the expectation of a statistic of the
j-th visited list (j ≤ m) is the j-th iterate of joinStep.
theorem exists_finset_randomJoining (m : ℕ) (b : List ℕ) :
∃ F : Finset (List (List ℕ)), (∀ ps ∉ F, randomJoining m b ps = 0) ∧
∀ j ≤ m, ∀ g : List ℕ → ℝ,
∑ ps ∈ F, randomJoining m b ps * g (ps.getD j []) = (joinStep^[j] g) b
noncomputable def randomJoining : ℕ → List ℕ → List (List ℕ) → ℝ
| 0, b, ps => if ps = [b] then 1 else 0
| m + 1, b, ps =>
if (blockPairs b.length).card = 0 then (if ps = [b] then 1 else 0)
else match ps with
| [] => 0
| b' :: rest =>
if b' = b then
(∑ p ∈ blockPairs b.length, randomJoining m (mergeBlocks b p.1 p.2) rest) /
(blockPairs b.length).card
else 0
MathlibExt/Biology/PopulationGenetics/Coalescent/SiteFrequencySpectrum.lean:84
noncomputable def joinStep (g : List ℕ → ℝ) (b : List ℕ) : ℝ :=
if (blockPairs b.length).card = 0 then g []
else (∑ p ∈ blockPairs b.length, g (mergeBlocks b p.1 p.2)) / (blockPairs b.length).card
MathlibExt/Biology/PopulationGenetics/Coalescent/SiteFrequencySpectrum.lean:100
Uses 4 · Used by 2
Order-free statistics and uniform compositions
TheoremExpected number of blocks of size i#✓fully proved
Expected number of blocks of size i when k blocks remain, for random joining
started from n singletons: k · C(n-i-1, k-2) / C(n-1, k-1) (Fu 1995; Griffiths–Tavaré
1998). Equivalently, a given one of the k blocks has size i with probability
C(n-i-1, k-2) / C(n-1, k-1).
theorem joinStep_iterate_count {n k i : ℕ} (hk : 2 ≤ k) (hkn : k ≤ n) (hi : 1 ≤ i)
(hin : i < n) :
(joinStep^[n - k] (fun b => (b.count i : ℝ))) (List.replicate n 1) =
k * ((n - i - 1).choose (k - 2) : ℝ) / ((n - 1).choose (k - 1) : ℝ)
noncomputable def joinStep (g : List ℕ → ℝ) (b : List ℕ) : ℝ :=
if (blockPairs b.length).card = 0 then g []
else (∑ p ∈ blockPairs b.length, g (mergeBlocks b p.1 p.2)) / (blockPairs b.length).card
MathlibExt/Biology/PopulationGenetics/Coalescent/SiteFrequencySpectrum.lean:100
Uses 3 · Used by 2
The expected branch length E L_i = 2 / i
TheoremExpected block counts weighted by coalescence times sum to \(2/i\)#✓fully proved
Summing the expected block counts against E T_k = (k.choose 2)⁻¹ gives 2 / i:
∑_{k=2}^{n} k C(n-i-1, k-2) / C(n-1, k-1) · (k.choose 2)⁻¹ = 2 / i for 1 ≤ i < n.
theorem sum_blockCount_mul_inv_choose_two {n i : ℕ} (hi : 1 ≤ i) (hin : i < n) :
∑ k ∈ Icc 2 n, k * ((n - i - 1).choose (k - 2) : ℝ) / ((n - 1).choose (k - 1) : ℝ) *
((k.choose 2 : ℕ) : ℝ)⁻¹ = 2 / i
The expected site frequency spectrum
TheoremExpected site frequency spectrum#✓fully proved
Expected site frequency spectrum (Fu 1995). Let the genealogy of a sample of n gene
copies consist of Kingman holding times K.T and an independent block-size history τ with
the random-joining law. If, given the genealogy (the σ-algebra generated by τ and the
holding times), the number ξ i of sites whose derived allele is carried by exactly i
copies is integrable and has conditional mean (θ / 2) L_i (neutral infinite-sites mutation,
θ = 4Nu), then E ξ_i = θ / i for 1 ≤ i ≤ n - 1. (Integrability matters: Mathlib's
conditional expectation of a non-integrable function is 0.)
theorem integral_sfs {μ : Measure Ω} [IsProbabilityMeasure μ] {n : ℕ}
(K : KingmanTimes μ n) (hT : ∀ k, Measurable (K.T k))
{τ : Ω → List (List ℕ)} (hτ : MeasurableSpace.comap τ ⊤ ≤ mΩ)
(hτ_law : ∀ ps, μ.real (τ ⁻¹' {ps}) =
randomJoining (n - 1) (List.replicate n 1) ps)
(h_indep : Indep (MeasurableSpace.comap τ ⊤)
(MeasurableSpace.comap (fun ω k => K.T k ω) MeasurableSpace.pi) μ)
{θ : ℝ} {ξ : ℕ → Ω → ℝ} (hξ_int : ∀ i, Integrable (ξ i) μ)
(hξ : ∀ i, μ[ξ i | MeasurableSpace.comap τ ⊤ ⊔
MeasurableSpace.comap (fun ω k => K.T k ω) MeasurableSpace.pi] =ᵐ[μ]
fun ω => θ / 2 * sfsLength n i (τ ω) (fun k => K.T k ω))
{i : ℕ} (hi : 1 ≤ i) (hin : i < n) :
μ[ξ i] = θ / i
structure KingmanTimes (μ : Measure Ω) (n : ℕ) where
/-- `T k` is the holding time while the genealogy has `k` lineages. -/
T : ℕ → Ω → ℝ
/-- `T k` is exponential with rate `k.choose 2`, the number of pairs of lineages. -/
hasLaw : ∀ k ∈ Icc 2 n, HasLaw (T k) (expMeasure (k.choose 2 : ℕ)) μ
/-- The holding times are mutually independent. -/
indep : iIndepFun (fun k : Icc 2 n => T k) μ
MathlibExt/Biology/PopulationGenetics/Coalescent/Kingman.lean:115
noncomputable def randomJoining : ℕ → List ℕ → List (List ℕ) → ℝ
| 0, b, ps => if ps = [b] then 1 else 0
| m + 1, b, ps =>
if (blockPairs b.length).card = 0 then (if ps = [b] then 1 else 0)
else match ps with
| [] => 0
| b' :: rest =>
if b' = b then
(∑ p ∈ blockPairs b.length, randomJoining m (mergeBlocks b p.1 p.2) rest) /
(blockPairs b.length).card
else 0
MathlibExt/Biology/PopulationGenetics/Coalescent/SiteFrequencySpectrum.lean:84
def sfsLength (n i : ℕ) (ps : List (List ℕ)) (T : ℕ → ℝ) : ℝ :=
∑ k ∈ Icc 2 n, ((ps.getD (n - k) []).count i : ℝ) * T k
MathlibExt/Biology/PopulationGenetics/Coalescent/SiteFrequencySpectrum.lean:107
Sources
- Fu (1995); Griffiths & Tavaré (1998). Proves the deferred entry
integral_sfs.
Uses 9 · Used by 1
Loading dependency explorer…
The dependency explorer needs JavaScript and its graph data.