Skip to content

Expected site frequency spectrum

9 fully proved

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].

def mergeBlocks (b : List ℕ) (i j : ℕ) : List ℕ :=
  ((b.eraseIdx j).eraseIdx i) ++ [b.getD i 0 + b.getD j 0]
Used by 4

DefinitionUnordered pairs of blocks#✓fully proved

The unordered pairs of positions i < j among k blocks.

def blockPairs (k : ℕ) : Finset (ℕ × ℕ) :=
  (range k ×ˢ range k).filter fun p => p.1 < p.2
Used by 4

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
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
Uses 2 · Used by 4

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.

def sfsLength (n i : ℕ) (ps : List (List ℕ)) (T : ℕ → ℝ) : ℝ :=
  ∑ k ∈ Icc 2 n, ((ps.getD (n - k) []).count i : ℝ) * T k
Used by 3

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
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) : ℝ)
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
Used by 2

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
Sources
  • Fu (1995); Griffiths & Tavaré (1998). Proves the deferred entry integral_sfs.
Uses 9 · Used by 1

The dependency explorer needs JavaScript and its graph data.

Download the complete graph data.