@swamidass/# TMR4A

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