Skip to content

Skyline plot estimates of population size

12 fully proved

Coalescent

StructureCoalescent intervals under the classic skyline model#✓fully proved

The coalescent intervals of a sample of n genes under the classic skyline model with piecewise-constant population-size parameter θ k while there are k lineages: T k is exponential with rate C(k, 2) / θ k, independently over 2 ≤ k ≤ n.

structure SkylineTimes (μ : Measure Ω) (n : ℕ) (θ : ℕ → ℝ) where
  /-- `T k` is the length of the interval with `k` lineages. -/
  T : ℕ → Ω → ℝ
  /-- The population-size parameters are positive. -/
  θ_pos : ∀ k ∈ Icc 2 n, 0 < θ k
  /-- `T k` is exponential with rate `C(k, 2) / θ k`. -/
  hasLaw : ∀ k ∈ Icc 2 n, HasLaw (T k) (expMeasure (((k.choose 2 : ℕ) : ℝ) / θ k)) μ
  /-- The interval lengths are mutually independent. -/
  indep : iIndepFun (fun k : Icc 2 n => T k) μ
Used by 4

DefinitionLog-likelihood of one coalescent interval in the skyline model#✓fully proved

The log-likelihood of θ for one coalescent interval of length w with c = C(k, 2) pairs of lineages: the log-density log (c / θ) - c w / θ of Exp(c / θ) at w.

def skylineLogLik (c w θ : ℝ) : ℝ := Real.log (c / θ) - c * w / θ
Used by 4

TheoremThe skyline estimate maximizes the interval likelihood#✓fully proved

The skyline estimate maximizes the interval likelihood (Pybus, Rambaut & Harvey 2000): for c, w, θ > 0, ℓ(θ) ≤ ℓ(c w).

theorem skylineLogLik_le {c w θ : ℝ} (hc : 0 < c) (hw : 0 < w) (hθ : 0 < θ) :
    skylineLogLik c w θ ≤ skylineLogLik c w (c * w)
Uses 1 · Used by 2

TheoremUniqueness of the skyline maximizer#✓fully proved

For c, w, θ > 0, ℓ(θ) = ℓ(c w) iff θ = c w.

theorem skylineLogLik_eq_iff {c w θ : ℝ} (hc : 0 < c) (hw : 0 < w) (hθ : 0 < θ) :
    skylineLogLik c w θ = skylineLogLik c w (c * w) ↔ θ = c * w
Uses 1 · Used by 1

TheoremThe skyline estimates jointly maximize the likelihood#✓fully proved

For positive c i, w i, θ i over a finite set of intervals, ∑ ℓ_i(θ_i) ≤ ∑ ℓ_i(c_i w_i).

theorem sum_skylineLogLik_le {ι : Type*} (s : Finset ι) {c w θ : ι → ℝ}
    (hc : ∀ i ∈ s, 0 < c i) (hw : ∀ i ∈ s, 0 < w i) (hθ : ∀ i ∈ s, 0 < θ i) :
    ∑ i ∈ s, skylineLogLik (c i) (w i) (θ i) ≤ ∑ i ∈ s, skylineLogLik (c i) (w i) (c i * w i)
Uses 2

TheoremUniqueness of the joint maximizer#✓fully proved

Equality in sum_skylineLogLik_le holds iff θ i = c i w i for every interval.

theorem sum_skylineLogLik_eq_iff {ι : Type*} (s : Finset ι) {c w θ : ι → ℝ}
    (hc : ∀ i ∈ s, 0 < c i) (hw : ∀ i ∈ s, 0 < w i) (hθ : ∀ i ∈ s, 0 < θ i) :
    ∑ i ∈ s, skylineLogLik (c i) (w i) (θ i) =
        ∑ i ∈ s, skylineLogLik (c i) (w i) (c i * w i) ↔ ∀ i ∈ s, θ i = c i * w i
Uses 3

DefinitionClassic skyline estimate of population size#✓fully proved

The classic skyline estimate θ̂_k = C(k, 2) w_k for the interval with k lineages.

def skylineEstimate (k : ℕ) (w : ℝ) : ℝ := ((k.choose 2 : ℕ) : ℝ) * w
Used by 3

Skyline Times

TheoremThe skyline estimate is unbiased#✓fully proved

E[C(k, 2) T k] = θ k.

theorem integral_skylineEstimate {k : ℕ} (hk : k ∈ Icc 2 n) :
    μ[fun ω => skylineEstimate k (S.T k ω)] = θ k
Uses 4

TheoremVariance of the skyline estimate#✓fully proved

Var[C(k, 2) T k] = θ k².

theorem variance_skylineEstimate {k : ℕ} (hk : k ∈ Icc 2 n) :
    Var[fun ω => skylineEstimate k (S.T k ω); μ] = θ k ^ 2
Uses 4

TheoremThe skyline estimates are independent#✓fully proved

The skyline estimates are independent across intervals.

theorem iIndepFun_skylineEstimate :
    iIndepFun (fun k : Icc 2 n => fun ω => skylineEstimate k (S.T k ω)) μ
Uses 2

Coalescent

DefinitionCanonical probability space for the skyline model#✓fully proved

The canonical probability space for the skyline model: the product of the exponential laws Exp(C(k, 2) / θ k), 2 ≤ k ≤ n.

def skylinePiMeasure (n : ℕ) (θ : ℕ → ℝ) (hθ : ∀ k ∈ Icc 2 n, 0 < θ k) :
    Measure (Icc 2 n → ℝ) :=
  haveI : ∀ k : Icc 2 n, IsProbabilityMeasure (expMeasure (((k : ℕ).choose 2 : ℕ) / θ k)) :=
    fun k => isProbabilityMeasure_expMeasure
      (div_pos (KingmanTimes.choose_two_pos k.2) (hθ k k.2))
  Measure.pi fun k : Icc 2 n => expMeasure (((k : ℕ).choose 2 : ℕ) / θ k)
Used by 1

DefinitionCoordinate projections form a skyline model#✓fully proved

The coordinate projections of skylinePiMeasure n θ form a skyline model whenever the θ k are positive; in particular SkylineTimes is inhabited.

MathlibExt.Biology.PopulationGenetics.Coalescent.SkylineTimes.ofPi · 24 lines
def SkylineTimes.ofPi (n : ℕ) {θ : ℕ → ℝ} (hθ : ∀ k ∈ Icc 2 n, 0 < θ k) :
    SkylineTimes (skylinePiMeasure n θ hθ) n θ :=
  haveI : ∀ k : Icc 2 n, IsProbabilityMeasure (expMeasure (((k : ℕ).choose 2 : ℕ) / θ k)) :=
    fun k => isProbabilityMeasure_expMeasure
      (div_pos (KingmanTimes.choose_two_pos k.2) (hθ k k.2))
  { T := fun k ω => if h : k ∈ Icc 2 n then ω ⟨k, h⟩ else 0
    θ_pos := hθ
    hasLaw := fun k hk => by
      have hT : (fun ω : Icc 2 n → ℝ => if h : k ∈ Icc 2 n then ω ⟨k, h⟩ else 0)
          = fun ω => ω ⟨k, hk⟩ := by
        funext ω
        simp [hk]
      rw [hT]
      exact (measurePreserving_eval
        (fun i : Icc 2 n => expMeasure ((((i : ℕ).choose 2 : ℕ) : ℝ) / θ i)) ⟨k, hk⟩).hasLaw
    indep := by
      have hT : (fun k : Icc 2 n => fun ω : Icc 2 n → ℝ =>
          if h : (k : ℕ) ∈ Icc 2 n then ω ⟨k, h⟩ else 0) = fun k ω => id (ω k) := by
        funext k ω
        simp [k.2]
      rw [hT]
      exact iIndepFun_pi
        (μ := fun i : Icc 2 n => expMeasure ((((i : ℕ).choose 2 : ℕ) : ℝ) / θ i))
        (fun _ => aemeasurable_id) }
Uses 3

The dependency explorer needs JavaScript and its graph data.

Download the complete graph data.