repl.it
Python 2.7

In population genetics, what is the ratio between TMRCA and the time to most recent four ancestors (TMR4A)?

fork
loading
Files
  • main.py
main.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import numpy as np

def n_kingman_waits(n):
  """the expected wait times / Ne of the 1 through n coalescents. 
  The n-th element is the (n+1) coalecent time, and is the expected time / Ne
  for n+2 alleles to coalesce to n+1. Ne is the effective population size (alleles).
  """
  k = np.arange(n) + 1
  return 4. / (k * (k+1)) 

def n_kingman_ages(n):
  # get the expected ages / Ne
  return np.cumsum(n_kingman_waits(n)[::-1])[::-1]

# Any large number will do.  The coalescents ages converge quickly as n grows.
expected_age = n_kingman_ages(107) 

# What do we multiply TMRCA by to estimate TMR4A?
# 4th coalescent corresponds with TMR4A
print expected_age[3] / expected_age[0]
#>>> 0.242990654206

# What do we multiply the first 5 coalescents by to estimate TMR4A?
print expected_age[3] / expected_age[:6]
#>>> [ 0.24299065  0.49056604  0.74285714  1.          1.26213592  1.52941176]