# What proportion of the people who lived 1000 years ago have genetic descendants alive today?

For context, I've been wondering about this for a paper I'm writing (in philosophy). Really, I want to figure out the chances that someone alive today will end up still having descendants 1000 years in the future. But people from 1000 years ago should be a good approximation.

This is not a full answer, but I post it as one because the reasoning is too long to post it as a comment.

You can make a rough estimation using a few assumptions.

1) Let's say 1 generation = 25 years, so 1000 years are 40 generations.

2) Let's say that half of the population is able to leave offspring, so the probability of one individual having descendants = 0.5

3) Supposing that the average number of siblings is 3.

• The probability of not having descendants in the 2nd generation = 0.5

• The probability of not having descendants in the 3rd generation is 0.5*(1-0.5^3) = 0.4375

After the 3rd generation, it becomes a combinatorial problem that is too hard for me. I can tell that the probability of not having offsprings, giving that it already has, is going to 0 because of the number of offspring > 2

If this approach appears useful, I suggest you ask for a more complete solution on https://math.stackexchange.com/

• Thanks, this is great! – HW. May 23 '19 at 17:53

Expanding on @heracho's answer and Wikipedia, assume pi denotes the probability of having (exactly) i children, and that dm denotes the probability of extinction by the mth generation (note that this includes extinctions by all generations n with n < m). Then dm can be expressed as:

$$\displaystyle d_m = p_0 + p_1 d_{m-1} + p_2 d_{m-1}^2 + p_3 d_{m-1}^3 + \dots = \sum_{i=0}^\infty p_i d_{m-1}^i$$

with d0 = 0.

In other words, the probability of going extinct by the mth generation is equal to the probability of i children going extinct in m – 1 generations, summed over i. (The child lineages are assumed independent.) This recurrence relation allows you to calculate the extinction (and, by extension, survival) probability at any desired number of generations, if the values of pi are known. In simple cases it is also possible to work out the asymptotic value of dm by solving the fixed point equation dm = dm – 1 (see the Wikipedia page for an example).

In relation to your specific application, we can assume roughly 40 generations and use some empirical estimates of pi. The parity distribution of Denmark is given by Frejka (2008) as p0 = 0.127, p1 = 0.179, p2 = 0.443, p3 = 0.251 (where I simplify by assuming that people who have strictly more than 3 children actually have only 3 children). Plugging these parameters into the recurrence relation and iterating for 40 generations (see R code below), we get d40 = 0.17. Thus some 83% of the population can be assumed to have descendants in 1000 years.

I have taken Denmark as an example simply because it is the first country mentioned in Table 2 of Frejka (2008). However, to illustrate the dependence of the asymptotic value of dm on the parity parameters consider the case of Italy. From the same publication, we find p0 = 0.240, p1 = 0.235, p2 = 0.369, p3 = 0.156. With these values, we compute d40 = 0.41, obviously a very different situation.

I am not a mathematician, nor a biologist (nor a demographer), and the above logic obviously employs many ridiculous simplifications. It would be very interesting to know how all this pans out in a more realistic model.

# Reference

Frejka, Tomas (2008) Parity distribution and completed family size in Europe: Incipient decline of the two-child family model? Demographic Research, 19(4): 47–72. https://www.demographic-research.org/volumes/vol19/4/19-4.pdf

# R code

# get extinction probability (d) in generation m, given extinction
# probability in generation m-1, conditioned on probabilities
# of having 0, 1, 2, ... children (vector p)
extinctionprob <- function(d, p) {
d <- d^(0:(length(p)-1))
sum(d*p)
}

# iterate for 'ngen' generations
iterate <- function(ngen, p) {
d <- 0
for (i in 1:ngen) d <- extinctionprob(d, p)
d
}

# usage example:
# iterate(ngen=40, p=c(0.240, 0.235, 0.369, 0.156))