Skyline plot estimates of population size
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) μ
T : ℕ → Ω → ℝ
MathlibExt/Biology/PopulationGenetics/Coalescent/Skyline.lean:64
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 / θ
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)
def skylineLogLik (c w θ : ℝ) : ℝ := Real.log (c / θ) - c * w / θ
MathlibExt/Biology/PopulationGenetics/Coalescent/Skyline.lean:74
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
def skylineLogLik (c w θ : ℝ) : ℝ := Real.log (c / θ) - c * w / θ
MathlibExt/Biology/PopulationGenetics/Coalescent/Skyline.lean:74
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)
def skylineLogLik (c w θ : ℝ) : ℝ := Real.log (c / θ) - c * w / θ
MathlibExt/Biology/PopulationGenetics/Coalescent/Skyline.lean:74
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
def skylineLogLik (c w θ : ℝ) : ℝ := Real.log (c / θ) - c * w / θ
MathlibExt/Biology/PopulationGenetics/Coalescent/Skyline.lean:74
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
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
def skylineEstimate (k : ℕ) (w : ℝ) : ℝ := ((k.choose 2 : ℕ) : ℝ) * w
MathlibExt/Biology/PopulationGenetics/Coalescent/Skyline.lean:127
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
def skylineEstimate (k : ℕ) (w : ℝ) : ℝ := ((k.choose 2 : ℕ) : ℝ) * w
MathlibExt/Biology/PopulationGenetics/Coalescent/Skyline.lean:127
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 ω)) μ
def skylineEstimate (k : ℕ) (w : ℝ) : ℝ := ((k.choose 2 : ℕ) : ℝ) * w
MathlibExt/Biology/PopulationGenetics/Coalescent/Skyline.lean:127
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)
lemma choose_two_pos {k : ℕ} (hk : k ∈ Icc 2 n) : (0 : ℝ) < ((k.choose 2 : ℕ) : ℝ)
MathlibExt/Biology/PopulationGenetics/Coalescent/Kingman.lean:134
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) }
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) μ
MathlibExt/Biology/PopulationGenetics/Coalescent/Skyline.lean:62
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)
MathlibExt/Biology/PopulationGenetics/Coalescent/Skyline.lean:163
lemma choose_two_pos {k : ℕ} (hk : k ∈ Icc 2 n) : (0 : ℝ) < ((k.choose 2 : ℕ) : ℝ)
MathlibExt/Biology/PopulationGenetics/Coalescent/Kingman.lean:134
Loading dependency explorer…
The dependency explorer needs JavaScript and its graph data.