Modelling the age distribution of longevity leaders
Csaba Kiss László Németh Bálint Vető


 Real-time diagnostics  on a  QKD link via
QBER Time Series Analysis 
     Dimitris Syvridis
    MOX - Dipartimento di Matematica “F. Brioschi”, Politecnico di
Milano, via Bonardi 9, 20133 Milan, Italy

    ^1 arashandrea.roknian@mail.polimi.it

    ^2 anna.scotti@polimi.it

    ^3 alessio.fumagalli@polimi.it
    September 9, 2024
===============================================================================================================================================================================================================================================





empty



§ INTRODUCTION


Exceptionally long human lifespans are one of the cornerstones of demography and mortality research.
Studying the group of record holders may reveal not only the underlying mortality mechanism of a population but also potentially shed some light on the future developments of human longevity.
As life expectancy increases <cit.> with deaths shifting to older ages, the distribution of deaths at the oldest-old ages <cit.> gains more interest of demographers, actuaries and decision makers of numerous disciplines.

The pattern of adult human mortality has already been described <cit.> but there is still a debate about the exact distribution of deaths at adult ages. More details on the necessity of a correctly specified model for the underlying mortality process and its impact on further research are discussed in <cit.>. Numerous publications review the ongoing discussion on the existence of a mortality plateau <cit.>, and the levelling-off of adult human death rates at the oldest ages is supported by the findings in <cit.> while others cast some doubt on this observation <cit.>. If there is a mortality plateau then the distribution of deaths at the oldest-old ages must be gamma-Gompertz and human lifespan can increase further without any maximum <cit.>. Further studies discuss the existence of a limit to human lifespan with more focus on the extreme value distribution aspect of the deaths at the oldest-old ages for various populations and different lifetime distributions <cit.>. These models can be helpful in determining the plausibility of longevity leaders as well.

We contribute to this discussion by proposing a stochastic model to describe the evolution of the age of the world's oldest person. Based on our estimates the model provides a good fit to the titleholder data since 1955, collected by the Gerontology Research Group <cit.>.
With the model results, it is possible to predict the age of the oldest person in the world in the future.
When should we expect to see the next Jeanne Calment, the supercentenarian with the longest human lifespan ever documented?
Will her record ever be surpassed?
Our results provide a prediction for the age distribution of the record holder in the coming decades to answer these questions.




§ RESULTS


Our model describes the evolution of the age of the oldest living person under the following assumptions.
We assume that the births of individuals follow a Poisson process with time-dependent intensity <cit.>.
The lifespans of individuals in the population are independent and their distribution may depend on the date of birth.
Then the age of the record holder in the population evolves in time as a Markov process with explicit transition probabilities.
As the first main result of this paper, we explicitly compute the distribution of the age of the record holder for any given birth rate parameter and lifespan distribution.
The detailed mathematical description and the properties of the general model are described in Section <ref>.

We apply our general result to the case which approximates the human birth rate over the world and the human lifespan.
We specify the intensity function of the Poisson process of births to have an exponential growth in time.
The underlying force of mortality is chosen so that it follows an extension of the Gompertz mortality model <cit.>
and the lifespan distribution of individuals is given by the gamma–Gompertz () distribution with time-dependent parameters.
This distribution adequately captures the slowing down of senescence mortality at the oldest old ages.
Given the growth parameters of the birth rate, we fit the model parameters to the statistics of the oldest person titleholder data using maximum likelihood method.
The optimal parameters of the model fit well to the data.
It shows in particular that the age of the oldest person alive increases over time, and it will most likely increase further in the future.
We compute the expected value and a confidence interval for the age of the world's oldest person using the fitted model parameters for each year between 1955 and 2019 shown by the green curves on Figure <ref>.
The detailed discussion of the model specification, likelihood calculations as well as the parameter fitting are given in Section <ref>.
Section <ref> contains calculations related to the gamma–Gompertz–Makeham generalization of the gamma–Gompertz distribution.




Our results enable us to predict the age distribution of the world's oldest person at future time points.
We compute the probability density of the age of the world's oldest person in different years not only in the past but also in the future.
These densities are shown in Figure <ref>.
When comparing the age distribution of the oldest person in the world in different years to the age of Jeanne Calment at her death, we find that on Jan 1st 2060 we can expect that the age of the world's oldest person will exceed her age with probability around 0.5.
This also means that with high probability her age record will already be broken by that time.



In Figure <ref> two extreme outliers with unexpectedly long lifetimes can be observed. Jeanne Calment died at the age of 122.45 years in 1997, and Sarah Knauss died at the age 119.27 in 2000.
In our model, the probability of observing an age greater than or equal to their actual age at the time of their death is 0.000286 for Calment and 0.0116 for Knauss.
See details in Subsection <ref>.

The fact that Calment and Knauss are outliers among the oldest old in the world became even more evident when we performed a backtesting of our model. We estimated the parameters based on the data on the world's oldest person between 1955 and 1988 where the ending date is the time when Calment became the world's oldest person. The model-estimated mean and confidence interval of the world's oldest person using the full data and the partial dataset (before Calment) are shown in Figure <ref> by green and red, respectively.
The estimate using the data until 1988 is less reliable after 2000 which is shown by the fact that the observed data is out of the confidence interval in the majority of the time after 2000.
When we compare the two confidence intervals, we can conclude that, based on the data before 1988, Calment and Knauss already had extremely high ages at their deaths.
Adding the remaining data set, the estimated mean age of the world's oldest person becomes lower.
Hence we can conclude that now we can consider the ages of the two outliers between 1988 and 2000 at their death to be more extreme than based on the information available in 1988.

The other important observable is the reign length of a record holder.
The numerical value of the expected reign length with our estimated model parameters is 1.195 in 1955 and it is 1.188 in 2019.
The empirical value of the reign length is 1.008 which is not much less than the model-based estimate.

Our approach to studying the age of the oldest old is completely novel because it takes into account jointly the age and the time of birth of individuals.
Although the age and the reign length of the world's oldest person depend in a complex non-linear way on the total lifespan and the time of birth of supercentenarians, we compute explicitly the probability distribution of the age of the oldest person.
Hence, the performance of our predictions cannot be directly compared to previous results in the literature because in the usual approach, the oldest person in each cohort is considered separately, and it is not relevant whether this person was ever the oldest in the population, see e.g. the extreme value method in <cit.>.
Our model contributes to the mathematical understanding of the evolution of the oldest individual, which is the extra benefit compared to a prediction using the trend in the data, e.g. in linear regression.
In this way, we not only observe but also prove mathematically that the dynamics of the birth process and that of the lifespan distribution which we consider in this paper necessarily imply the increase of the expected age of the world's oldest person.



§ MATHEMATICAL MODEL FOR THE AGE OF THE OLDEST PERSON



In this section, we provide the mathematical definition of a general model for the age of the world's oldest person, where the births of individuals follow a Poisson process and their lifespans are independent.
Under the assumptions of the general model, the age of the record holder in the population evolves in time as a Markov process with explicit transition probabilities.
In Subsections <ref>–<ref>, we describe the exact distribution of the age of the record holder in this generality
for any given birth rate parameter and lifespan distribution using the two-dimensional representation of the age process of the oldest person.
In the time-homogeneous case with constant birth rate and identical lifespan distributions the reign length distribution of a record holder is computed in Subsections <ref>–<ref>.
We explain the role of the entry age parameter in Subsection <ref>.



 §.§ Model description and two-dimensional representation



The model is formally defined as follows.
Let λ(t) be the birth rate parameter which depends on time and let F_t and f_t be a family of cumulative distribution functions and density functions corresponding to non-negative random variables which are also time-dependent.
We assume that individuals are born according to a Poisson point process at rate λ(t) and that the lifespan of an individual born at time t is given by F_t so that lifespans are independent for different individuals.

Let Y_t denote the age of the oldest person in the population at time t.
The process (Y_t:t∈) is Markovian.
The Markov property holds because at any time t the history of the process (Y_s:s≤ t) provides information about the lifetime of individuals born before the current record holder
while any transition of (Y_s:s≥ t) depends only on the lifetime of the current record holder and of those born after them.

The evolution of the Markov process Y_t is the following.
It has a deterministic linear growth with slope 1 due to the ageing of the current record holder.
This happens until the death of the record holder.
Additionally, given that Y_t-=lim_s↑ tY_s=y for some t with y>0,
the process has a downward jump at time t at rate f_t(y)/(1-F_t(y)) which is the hazard rate of the distribution F_t at y.
This corresponds to the possibility that the record holder dies at time t which happens at rate f_t(y)/(1-F_t(y)).
The conditional distribution of the jump is given by

    (Y_t<x | Y_t-=y,Y_t<y)
    =exp(-∫_x^yλ(t-u)(1-F_t-u(u)) ụ)

for all x>0.
The jump distribution in (<ref>) has an absolutely continuous part supported on [0,y] with density

    j_y,t(x)=exp(-∫_x^yλ(t-u)(1-F_t-u(u)) ụ)
    λ(t-x)(1-F_t-x(x))

and a point mass at 0 with probability

    a_y,t=exp(-∫_0^yλ(t-u)(1-F_t-u(u)) ụ).

As we shall see in the relevant parameter regime the probability a_y,t of the point mass at 0 is negligible.
The transition formula in (<ref>) can be proven using the description below.

We introduce a two-dimensional representation of the process Y_t as follows.
Let Λ={(t_i,x_i):i∈ I} be a marked Poisson process in ×_+
where {t_i:i∈ I} forms a Poisson point process on  with intensity λ(t)
and x_i≥0 is sampled independently for each i∈ I according to the distribution F_t_i.
The point (t_i,x_i) represents an individual born at time t_i with lifespan x_i for all i∈ I, that is, the individual i is alive in the time interval [t_i,t_i+x_i) and their age at time t is t-t_i if t∈[t_i,t_i+x_i).
Hence the marked Poisson process Λ contains all relevant information about the age statistics of the population at any time.
In particular the age of the oldest person Y_t can be expressed in terms of Λ as

    Y_t=max{(t-t_i)𝕀_{t∈[t_i,t_i+x_i)}:i∈ I}.

where the indicator 𝕀_{t∈[t_i,t_i+x_i)} is 1 exactly if the ith person is alive at time t.

The transition distribution formula (<ref>) can be seen using the two-dimensional representation as follows.
Given that the current record holder dies at time t at age y the event {Y_t<x} means that nobody with age between x and y can be alive at time t.
This event can be equivalently characterized in terms of the Poisson process of birth at rate λ(·) thinned by the probability that the person is still alive at time t.
Indeed the event {Y_t<x} can be expressed as a Poisson process of intensity at time t-u given by λ(t-u)(1-F_t-u(u)) for u∈[x,y] not having any point in the time interval [t-y,t-x].
This probability appears exactly on the right-hand side of (<ref>).
In other words, for any u∈[x,y] people are born at time t-u at rate λ(t-u).
On the other hand, the probability for a person born at time t-u to be alive at time t (that is, at age u) is 1-F_t-u(u).
See Figure <ref> for illustration.





 §.§ Exact distribution of the oldest person's age process



We assume that all birth events are already sampled on (-∞,t] together with the corresponding lifespans.
Then the distribution of Y_t can be computed explicitly for all t∈ using the two-dimensional representation.
For all t∈ the density

    h_t(x)=exp(-∫_x^∞λ(t-u)(1-F_t-u(u)) ụ)
    λ(t-x)(1-F_t-x(x))

for all x>0 and the point mass at 0

    m_t=exp(-∫_0^∞λ(t-u)(1-F_t-u(u)) ụ)

characterize the distribution of Y_t which can be seen as follows.
We mention that the point mass at 0 is negligible in the application.

Similarly to the proof of the transition formula in (<ref>) the event {Y_t<x} for any x>0 is the same as the event that nobody with age at least x is alive at time t.
We express this event in terms of the Poisson process of birth at rate λ(·) thinned by the probability that the person is still alive at time t.
The event {Y_t<x} means that a Poisson process of intensity at time t-u given by λ(t-u)(1-F_t-u(u)) for u≥ x does not have any point in (-∞,t-x] yielding

    (Y_t<x)=exp(-∫_x^∞λ(t-u)(1-F_t-u(u)) ụ).

In other words for any u≥ x individuals are born at time t-u at rate λ(t-u).
A person born at time t-u is alive at time t at age u with probability 1-F_t-u(u).
(<ref>)–(<ref>) follow by differentiation in (<ref>) and by taking the x→0 limit.
See Figure <ref> for illustration.





 §.§ Homogeneous model



The exact computation of the reign length distribution (see Subsection <ref>) can only be performed in a special case of our general model described in Subsection <ref>.
We introduce this special case as the homogeneous model where individuals are born at the times of a Poisson process of constant rate λ=λ(t)>0 for all t.
The lifespan of individuals are independent and identically distributed with a fixed density f=f_t and cumulative distribution function F=F_t for all t which does not depend on time.

In the homogeneous model the jump distribution given in (<ref>)–(<ref>) simplifies to

    j_y(x)   =j_y,t(x)=exp(-λ∫_x^y(1-F(u)) ụ)λ(1-F(x)),
    
    a_y   =a_y,t=exp(-λ∫_0^y(1-F(u)) ụ).

The distribution of Y_t does not depend on time in this case hence it is a stationary distribution as well.
The formulas for the density of Y_t and point mass at 0 reduce in the homogeneous case to

    h(x)   =exp(-λ∫_x^∞(1-F(u)) ụ)λ(1-F(x)),
    
    m   =exp(-λ∫_0^∞(1-F(u)) ụ)

where the integral ∫_0^∞(1-F(u)) ụ is equal to the expected lifespan.

The equilibrium condition for the homogeneous density h can be written as

    h'(x)+h(x)f(x)/1-F(x)-∫_x^∞ h(y)f(y)/1-F(y)j_y(x) ỵ=0.

After differentiation and using the fact that

    /x̣j_y(x)=j_y(x)(λ(1-F(x))-f(x)/1-F(x))

one can derive from (<ref>) the second order differential equation

    h”(x)+(2f(x)/1-F(x)-λ(1-F(x)))h'(x)
    +(2f(x)^2/(1-F(x))^2+f'(x)/1-F(x))h(x)=0.

The point mass m at 0 satisfies

    mλ=∫_0^∞ h(y)f(y)/1-F(y)j_y(0) ỵ.




 §.§ The peaks process



In the homogeneous model the sequence of peaks in Y_t forms a discrete time Markov chain.
By peak we mean a local maximum of Y_t with value being equal to the lifespan of the last record holder.
Each time the oldest person dies the process Y_t has a peak with a downward jump following it.
Let Z_n denote the age of record holders at which they die which are the values of the peaks of the process Y_t.
The sequence Z_n forms a discrete time Markov chain.
The Markov property follows by the fact that ages at death of previous record holders only give information on people born before the current record holder but transitions depend on the lifespan of the current record holder and that of people born after them.

The stationary density of Z_n is given by

    z(x)=f(x)exp(-λ∫_x^∞(1-F(u)) ụ)/∫_0^∞ f(y)exp(-λ∫_y^∞(1-F(u)) ụ)ỵ.

The formula can be seen as follows.
To have a record holder who dies at age x there has to be a person who has lifespan x which gives the factor f(x) in the numerator on the right-hand side of (<ref>).
The exponential factor is by the two-dimensional representation equal to the probability that no people born before the record holder who just died can be alive at the time the record holder dies.
The denominator on the right-hand side of (<ref>) makes z(x) a probability density function.

The density of Z_n can also be characterized by the following description.
It satisfies the integral equation

    z(x)=∫_0^∞ z(w)∫_0^min(x,y)j_w(y)f(x)/1-F(y) ỵ ẉ+∫_0^∞ z(w)a_wf(x) ẉ

which comes from the possible transitions of the peak process as follows.
If the previous record holder had a total lifetime w∈[0,∞) then at the death the process Y_t jumps down to some value y at rate j_w(y) or to 0 with probability a_w.
The density of the age at which a person dies who becomes a record holder at age y is f(x)/(1-F(y)).
From the integral equation in (<ref>) one can derive the second order differential equation for the function g(x)=z(s)/f(x) given by

    g”(x)-λ(1-F(x))g'(x)+λ f(x)g(x)=0

which is satisfied by g(x)=cexp(-λ∫_x^∞(1-F(u)) ụ) in accordance with (<ref>).



 §.§ Reign length distribution



In the homogeneous model, let W_n denote the reign length of the nth record holder, that is, the time length for which this person is the oldest person of the population.
The density of the random reign length is given by

    r(w)=∫_0^∞ h(y)f(y+w) ỵ+m∫_0^wf(z)λ e^-λ(w-z)f(z) ẓ/∫_0^∞ h(y)(1-F(y)) ỵ+m.

The density formula in (<ref>) can be derived based on the stationary density of the peaks process given by (<ref>) as follows.

It holds for the density of the reign length that

    r(w)=∫_0^∞ z(x)(∫_0^xj_x(y)f(y+w)/1-F(y) ỵ
    +a_x∫_0^wf(z)λ e^-λ(w-z) ẓ)x̣

based on the decomposition with respect to the previous value of the peaks process Z_n.
The integral ∫_0^wf(z)λ e^-λ(w-z) ẓ is the density of the convolution of the density f with an independent exponential distribution of parameter λ.
On the right-hand side of (<ref>), one can use the definitions of the density z given by (<ref>) and the homogeneous jump distribution j_x and a_x given by (<ref>).
Then in the numerator after the exchange of the order of integrations in the first term and by using the formula for the stationary distribution given by (<ref>) one gets the numerator of (<ref>).
In the denominator one can use the equality

    ∫_0^∞ f(y)exp(-λ∫_y^∞(1-F(u)) ụ)ỵ
    =∫_0^∞ h(y)(1-F(y)) ỵ+m.

which follows by integration by parts.

Note also that the density is not equal to the remaining reign length of Y_t under the stationary distribution because it would involve the integral ∫_0^∞ h(y)f(y+w)/(1-F(y)) ỵ in place of the first term in the numerator on the right-hand side of (<ref>).



 §.§ Entry age parameter



Next we introduce another parameter which we call the entry age and we denote it by E.
As opposed to our original model we consider individuals as being born at age E at a modified birth rate and with a modified lifespan distribution.
As a result we obtain a model to the age of the world's oldest person all values of the entry age parameter E≥0 and we can fit the model parameters with different values of the entry age.

For any value E≥0 of the entry age we denote by λ_E(t) the rate at which people reach the age E at time t, that is,

    λ_E(t)=λ(t-E)(1-F_t-E(E))

because the new birth process of rate λ_E(t) is obtained by an inhomogeneous thinning of the original Poisson process of the birth events.
The lifespan distribution of those born at time t with age E becomes the remaining lifetime distribution at age E.
The modified cumulative distribution function and density are given by

    F_t^E(x)=F_t-E(x+E)-F_t-E(E)/1-F_t-E(E),    
    f_t^E(x)=f_t-E(x+E)/1-F_t-E(E).




§ MODEL SPECIFICATION AND PARAMETER FITTING



In Subsection <ref> we specify the general model introduced and discussed in Section <ref>, that is,
we assume that the birth rate parameter increases exponentially in time and that the lifespan distribution is given by the gamma–Gompertz–Makeham distribution with time-dependent parameters.
We provide the details of the computation of the likelihood as a function of the model parameters in Subsection <ref>.
We show the way to maximize the likelihood and how the optimal parameters can be found using the Nelder–Mead method in Subsection <ref>.
With these values of the parameters, the age of the world's oldest person and the reign length of the record holder can be computed as described in Subsection <ref>.



 §.§ Model specification: birth rate parameter and lifespan distribution



For the rest of the paper we specifiy our general model described in Section <ref>
to the following choice of the birth rate parameter and of the lifespan distribution.
We choose the value of the entry age to be E=0,30,60 and we fit the model parameters for all three values of E separately.

First we specify the intensity function of the Poisson process of births with an exponential growth in time.
For any of the three values of the entry age E we assume that the birth rate at age E is given by

    λ_E(t)=C_Ee^κ_Et

where the numerical values of the parameters C_E and κ_E are obtained
by linear regression of the logarithmic data of newborns, people at age 30 and 60 published by the United Nations since 1950.
We extrapolate the linear regression backwards in time and we use the numerical values shown in Table <ref>.

We assume that the underlying force of mortality is chosen so that the lifespan distribution of individuals follows the gamma–Gompertz () distribution with cumulative distribution function and density

    F^_a,b,γ(x)   =1-(1+aγ/b(e^bx-1))^-1/γ,
    
    f^_a,b,γ(x)   =ae^bx(1+aγ/b(e^bx-1))^-1-1/γ

for x≥0 where a,b,γ are positive parameters.
We mention that the gamma–Gompertz–Makeham () distribution differs from the  distribution by the presence of a non-negative extrinsic mortality parameter c which appears as an additive term in the force of mortality.
See (<ref>) for the definition of the  distribution.
In our model, we exclude the extrinsic mortality for the following two reasons.
Since the extrinsic mortality becomes irrelevant at high ages and we aim to model the front-end of the death distribution at the oldest-old ages,
we do not expect to obtain a reliable estimate on the extrinsic mortality using the data about the world's oldest person.
On the other hand, as explained later in Subsection <ref>, the likelihood maximization provides unrealistic lifespan distributions even for the  model if one tries to optimize in all the parameters at the same time.
In order for the algorithm to result in a distribution close to the actual human lifespan distribution, the number of model parameters had to be decreased.

For our model we suppose that in the  lifespan distribution, parameters b=b_E, the rate of aging and γ=γ_E, the magnitude of heterogeneity
are constants over time and that they only depend on the value of the entry age parameter E.
The parameter a, the initial level of mortality at the entry age for individuals born at time t, depends on time given by the exponentially decreasing function

    a_E(t)=K_Ee^-α_E(t-2000)

where the exponent α_E and the constant K_E only depends on the entry age E.
The reason for subtracting 2000 in (<ref>) is only technical, the numerical values of the parameters do not become tiny with this definition.

In the model with entry age E, we assume that the birth rate λ_E(t) is given by (<ref>) and we fit the gamma–Gompetz distribution with parameters b_E,γ_E and a_E(t) given by (<ref>) for the modified distribution function F_t^E(x) and density f_t^E(x) in (<ref>).
This means that we search for the best fitting values of the parameters α_E,K_E,b_E,γ_E which results in an approximation of the remaining lifetime distribution at the age E.



 §.§ Likelihood calculations



The aim of the maximum likelihood method is to give an estimate to the parameters α_E, K_E, b_E and γ_E for E=0,30,60 by finding those values for which the likelihood of the full sample is the largest.
The sample is obtained from the the historical data on the world's oldest person available in <cit.>.
We transform this information into a list of triples of the form (t_i,y_i,z_i) for i=1,…,n
where t_i is the ith time in the sample when the oldest person dies at age y_i and the new record holder has age z_i at time t_i.
Then the data has to satisfy the consistency relation t_i-z_i=t_i+1-y_i+1 since the two sides express the date of birth of the same person.

In the model with entry age E, the likelihood of the ith data point (t_i,y_i,z_i) given the previous data point is equal to

    f_t_i-y_i+E^E(y_i-E)/1-F_t_i-y_i+E^E(z_i-1-E)j^E_y_i,t_i(z_i)
    
    =f_t_i-y_i+E^E(y_i-E)/1-F_t_i-y_i+E^E(z_i-1-E)exp(-∫_z_i^y_iλ_E(t_i-u+E)(1-F_t_i-u+E^E(u-E)) ụ)
    λ_E(t_i-z_i+E)(1-F_t_i-z_i+E^E(z_i-E))

for all i=2,3,…,n except for i=1 in which case the 1-F_y_1-t_1+E^E(z_0-E) factor in the denominator is missing.
In (<ref>) above we use the transition probabilities of the model with entry age E given by

    j^E_y,t(x)=exp(-∫_x^yλ_E(t-u+E)(1-F_t-u+E^E(u-E)) ụ)λ_E(t-x+E)(1-F_t-x+E^E(x-E))

as a generalization of (<ref>).
The explanation of the left-hand side of (<ref>) is that the person died at time t_i at age y_i had age E at time t_i-y_i+E.
The previous data point ensures that this person has already reached age z_i-1 hence we condition their lifetime distribution on this fact.
The transition probabilities in (<ref>) are obtained similarly to (<ref>) with the difference that a person at age u with u∈[x,y] at time t had age E at time t-u+E.

Note that when computing the likelihood of the full data by multiplying the right-hand side of (<ref>) for different values of i the consistency relation of the data implies that the factor 1-F_t_i-z_i+E(z_i-E) of the ith term cancels with the factor 1-F_t_i+1-y_i+1+E(z_i-E) coming from the (i+1)st term.
Hence the log-likelihood of the full sample is given by

    l(α,K,b,γ)
       =∑_i=1^n(log f_t_i-y_i+E^E(y_i-E)
    -∫_z_i^y_iλ_E(t_i-u+E)(1-F_t_i-u+E^E(u-E)) ụ+logλ_E(t_i-z_i+E))
    +log(1-F_t_n-z_n-E^E(z_n+E))
       =∑_i=1^n(log f^_Ke^-α(t_i-y_i+E-2000),b,γ(y_i-E)
    -∫_z_i^y_iCe^κ(t_i-u+E)(1-F^_Ke^-α(t_i-u+E-2000),b,γ(u-E))ụ)
           +log(1-F^_Ke^-α(t_n-z_n+E-2000),b,γ(z_n-E))
    +nlog C+∑_i=1^nκ(t_i-z_i+E).

where we suppress the dependence of the parameters α,K,b,γ on the entry age.
Note that the last two terms do not depend on the parameters α,K,b,γ hence we can omit these terms in the maximization of the log-likelihood.



 §.§ Likelihood maximization



We implemented the calculation of the log-likelihood function l(α,K,b,γ) given by (<ref>) in Python.
We used numerical integration to obtain the integrals on the right-hand side of (<ref>).
We mention that the general integral formula in (<ref>) could not be used because the parameter a of the gamma–Gompertz distribution in the integrand depends on the integration variable on the right-hand side of (<ref>).

In order to maximize the value of the log-likelihood function l(α,K,b,γ) we applied the Nelder–Mead method
<cit.> which is already implemented in Python.
We mention that initially we used the gamma–Gompertz–Makeham distribution as lifespan distribution, see (<ref>) for the definition, which contains the extra parameter c to be fitted but it turned out that the number of model parameters has to be reduced.
The behaviour of the optimization algorithm in the five parameters α,K,b,c,γ using the gamma–Gompretz–Makeham model was very similar the case of four parameters α,K,c,γ in the gamma–Gompertz model.
Running the optimization in the full set of parameters (α,K,b,c,γ in the gamma–Gompertz–Makeham model or α,K,b,γ in the gamma–Gompertz model), it turned out that after a few rounds the parameter K started to decrease dramatically and reached values below 10^-10.
The resulting lifespan distribution seemed very unrealistic with almost no mortality before the age of 100.
This happened for all values of the entry age E=0,30,60.

We explain this phenomenon by the fact that historical data about the oldest person in the world only gives information about the behaviour of the lifespan distribution between the ages 107 and 123.
The simple optimization in the four parameters α,K,b,γ simultaneously yields an excellent fit for the tail decay of the lifespan distribution with the historical data but the result may be very far from the actual human lifespan.
This would limit the practical relevance of our results.

The mathematical reason for the fact that the four-parameter optimization does not result in a satisfactory approximation to the human lifespan distribution is the following.
In these cases, the optimization procedure diverges to those regimes of the parameter space _+^4 where the corresponding gamma–Gompertz distribution is degenerate.
One can prevent reaching these unrealistic combinations of parameters by reducing the amount of freedom in the optimization.
Hence we specify some of the parameters a priori and we perform the optimization in the remaining ones so that it provides a good fit to the data on the age of the oldest old as well as a realistic lifespan distribution.

We believe that the most robust of the four parameters of the  model is b which is the exponent in the time dependence of the mortality rate.
By setting the rate of aging b=0.09 the algorithm gives the optimal triple α,K,γ with the best likelihood which is very stable under changing the initial values of these parameters.
The running time is also very short.

The Nelder–Mead algorithm, being a numerical maximization method, heavily relies on the tolerance parameter, which determines the minimal improvement required for the algorithm to continue running. If this parameter is set too high, the algorithm might stop before reaching the optimum. Conversely, if set too low, the algorithm might take excessively long to converge. To address this, we drew inspiration from dynamic learning rate algorithms used in neural network training and developed the following meta-algorithm.

First, we run the Nelder–Mead optimization. Based on the improvement from the starting point, we dynamically adjust the tolerance factor, similar to how learning rates are modified during neural network training. We then run the optimization again, recalibrating the tolerance factor based on the observed improvement, and repeat the process. This iterative adjustment allows us to get closer to the optimum, a hypothesis supported by our practical experience with this meta-algorithm.
Following this meta-algorithm, only a few calls of the Nelder–Mead method is enough to reach the optimum.
The Python codes for the likelihood calculations as well as the Nelder–Mead optimization implemented to this problem are available in <cit.>.

The numerical values of the resulting parameters for the three choices of the entry age are shown in Table <ref>.
The survival probability functions with the parameters given in Table <ref> for individuals born in 2000 corresponding to the entry age E=0,30,60 are shown on Figure <ref> as a function of the age.
We also computed the optimal values of the parameters α,K,γ for other values of the rate of aging b as a sensitivity analysis.
The resulting parameter values for the choices b=0.11, b=0.13 and b=0.15 are shown in Table <ref>.



We mention as an alternative approach that scaling the parameters could enhance the optimization process, but this requires prior knowledge of the range within which the parameters vary. This range could be determined through our iterative application of the Nelder–Mead algorithm.



 §.§ Computation of the oldest person's age and of the reign length



We observe that the model with the parameters in Table <ref> fits well to the titleholder data.
We focus on two statistics of the process in order to support this observation about the comparison:
the age of the world's oldest person and the reign length of the record holder.
In the case of both statistics exact formulas are only available for the homogeneous model introduced in Subsection <ref> where the birth rate is constant as well as the lifespan distribution does not depend on time.
Hence we apply an approximation where the error is negligible compared to the difference from the statistics computed using the data.

In the general model the distribution of the age of the world's oldest person at time t is given by the density h_t(x) in (<ref>) and by the point mass m_t at 0 in (<ref>).
For the numerical computations, we ignore the point mass m_t which is below the round-off error in the numerical results.
The difficulty in computing the mean age of the oldest person at time t is that parameter a of the gamma–Gompertz–Makeham distribution function F_t-u in the exponent of (<ref>) also depends on the integration variable u.

In our approximation we fix the value of the parameter a of the  distribution in h_t(x) in (<ref>) to a value which is equal to a_0(t-d) in (<ref>) with some delay d.
The delay d is chosen so that the mean age of the oldest person computed using a_0(t-d) as parameter a for all times in the  distribution function in (<ref>) is equal to the same value d.
For a given t, this value of d can be obtained as the fixed point of the contraction map

    d↦∫_0^∞ xe^-∫_x^∞λ(t-u)(1-F_t-d(u)) ụλ(t-x)(1-F_t-d(x)) x̣

which provides a reasonable approximation for the mean age of the world's oldest person.
For the comparison with the data and for the prediction, we use the model with entry age E=0.
Hence in (<ref>), the function λ is given in (<ref>) with E=0 and F_t-d is the  distribution function with parameters given by the E=0 values in Table <ref> and with a=a_0(t-d) in (<ref>).
This approximation is reasonable because the distribution of the age of the oldest person is highly concentrated.

The fixed point of the map in (<ref>) as the expected age of the world's oldest person can be found in a few steps of iterations.
We show the result on Figure <ref>.
We applied 10 iterations using the fitted model parameters for each year between 1955 and 2019.
By computing the standard deviation of the age of the oldest person as well we obtain the mean and a confidence interval for the age.

The predictions for the age distribution of the world's oldest person in the future shown on Figure <ref>.
We obtained them by using the exact age distribution formula in Subsection <ref> along with the numerical values of the parameters C,κ,α,K,γ given in Tables <ref> and <ref> for entry age E=0.

In our model, the distribution function of the age of the world's oldest person is given in (<ref>) where the  distribution function can be substituted with the estimated parameter values at any time.
In this way, the probability of observing an age greater than or equal to Calment's or Knauss' actual age at the time of their death can be computed exactly.
The numerical values are 0.000286 for Calment and it is 0.0116 for Knauss.

The backtesting mentioned in the Results section is performed as follows.
We estimated the best parameter values with entry age 0 based on the reduced data on the world's oldest person between 1955 and 1988 where the ending date is the time when Calment became the world's oldest person.
The resulting parameters α=0.01516,K=0.00002064,γ=0.08413 are numerically not very far from the optimal parameters in Table <ref> but the difference is more visible on Figure <ref>.
The figure shows the model-based mean age and confidence interval for the age of the world's oldest person computed using the full data as well as the data until 1988.

For the reign length of record holders, we again used the expected age at a given time obtained as the fixed point of the iteration in (<ref>).
The numerical value of the expected reign length obtained from the iteration is 1.195 in 1955 and it is 1.188 in 2019.
The empirical value of the reign length is 1.008 computed from the data by dividing the total length of the time interval between 1955 and 2019 by the number of record holders.



§ METHODS



In this section, we provide supplementary information related to the main result of this paper.
We perform explicit computations with the gamma–Gompertz–Makeham model
and we express the integral of the survival function in terms of a hypergeometric function.

The cumulative distribution function and the density of the gamma–Gompertz–Makeham () distribution are given by

    F^_a,b,c,γ(x)   =1-e^-cx/(1+aγ/b(e^bx-1))^1/γ,
    
    f^_a,b,c,γ(x)   =e^-cx/(1+aγ/b(e^bx-1))^1+1/γc(b-aγ)+a(b+cγ)e^bx/b

for x≥0 where a,b,c,γ are positive parameters.
The positivity of parameters implies the finiteness of all moments and, in particular, the convergence of the integral of the survival function
∫_x^∞(1-F^_a,b,c,γ(u)) ụ.
In the homogeneous model, the integral of the survival function appears in the density of the distribution of Y_t in (<ref>) and in the stationary density of the peaks process in (<ref>).
We show below that in the gamma–Gompertz–Makeham model the integral of the survival function can be computed explicitly and it is given by

    ∫_x^∞(1-F^_a,b,c,γ(u))ụ
    =( b/a γ)^1/γe^-(c+b/γ)x/b/γ+c_2F_1(1/γ,1/γ+c/b;1+1/γ+c/b;a γ-b/a γe^-bx)

where _2F_1(a,b;c;z) is the hypergeometric function.
See 15.1.1 in <cit.> for the definition and properties.

We prove (<ref>) based on the following integral representation 15.3.1 in <cit.> of the hypergeometric function

    _2F_1(a,b,c,z)=Γ(c)/Γ(b) Γ(c-b)∫_0^1 t^b-1(1-t)^c-b-1(1-tz)^-a ṭ

which holds whenever (c)>(b)>0.
First we prove an identity for complex parameters α,β,δ which satisfy (α+β)>0 and we compute

    ∫_x^∞e^-β u/(1+δ e^u)^α ụ   =1/δ^α∫_x^∞e^-(α+β)u/(1+e^-u/δ)^α ụ
       =e^-(α+β)x/δ^α∫_0^1 y^α+β-1(1+e^-xy/δ)^-α ỵ
       =e^-(α+β)x/(α+β)δ^α·_2F_1(α,α+β;1+α+β;-e^-x/δ)

where we applied a change of variables y=e^x-u in the second equality above and we applied the hypergeometric identity (<ref>) in the last equality with a=α, b=α+β, c=1+α+β, z=-e^-x/δ together with the observation that with these values of the parameters the prefactor of the integral on the right-hand side of (<ref>) simplifies to α+β.
Note that the condition (c)>(b)>0 for (<ref>) to hold is satisfied by our assumption (α+β)>0 which also makes the integrals in (<ref>) convergent.

Next we show (<ref>) using (<ref>) as follows.
We write

    ∫_x^∞(1-F^_a,b,c,γ(u)) ụ   =1/(1-aγ/b)^1/γ∫_x^∞e^-cu/(1+aγ/b-aγe^bu)^1/γ ụ
       =1/(1-aγ/b)^1/γb∫_bx^∞e^-cv/b/(1+aγ/b-aγe^v)^1/γ ṿ
       =1/(1-aγ/b)^1/γbe^-(1/γ+c/b)bx/(1/γ+c/b)(aγ/b-aγ)^1/γ_2F_1(1/γ,1/γ+c/b;1+1/γ+c/b;a γ-b/a γe^-bx)

where we applied the change of variables v=bu in the second equality above and we used (<ref>) with α=1/γ, β=c/b, δ=aγ/(b-aγ) and with x replaced by bx.
The right-hand side of (<ref>) simplifies to that of (<ref>).

10
rm
url<#>1urlprefixURL doiprefixDOI: 


OepVau02
authorOeppen, J. & authorVaupel, J. W.
journaltitleBroken limits to life expectancy.
Science volume296,
  pages1029–1031 (year2002).

Can10
authorCanudas-Romo, V.
journaltitleThree measures of longevity: Time
  trends and record values.
Demography volume47,
  pages299–312 (year2010).

Vauetal21
authorVaupel, J. W., authorVillavicencio, F. &
  authorBergeron-Boucher, M.-P.
journaltitleDemographic perspectives on the
  rise of longevity.
Proceedings of the National Academy of
  Sciences volume118, pagese2019536118
  (year2021).

Gom25
authorGompertz, B.
journaltitleOn the nature of the function
  expressive of the law of human mortality, and on a new mode of determining
  the value of life contingencies.
Philosophical Transactions of the Royal Society
  of London volume115, pages513–583
  (year1825).

NemMis18
authorNémeth, L. & authorMissov, T. I.
journaltitleAdequate life-expectancy
  reconstruction for adult human mortality data.
Plos one volume13,
  pagese0198485 (year2018).

Misetal16
authorMissov, T. I., authorNémeth, L. &
  authorDańko, M. J.
journaltitleHow much can we trust life tables?
  sensitivity of mortality measures to right-censoring treatment.
Palgrave Communications
  volume2, pages1–10 (year2016).

Maietal21
authorMaier, H., authorJeune, B. &
  authorVaupel, J. W.
titleExceptional lifespans
  (publisherSpringer Nature, year2021).

Danetal23
authorDang, L. H. K. et al.
journaltitleThe question of the human mortality
  plateau.
Demographic Research
  volume48, pages321–338
  (year2023).

VijLeb17
authorVijg, J. & authorLe Bourg, E.
journaltitleAging and the Inevitable Limit to
  Human Life Span.
Gerontology volume63,
  pages432–434 (year2017).

Alvetal21
authorAlvarez, J.-A., authorVillavicencio, F.,
  authorStrozza, C. & authorCamarda, C. G.
journaltitleRegularities in human mortality
  after age 105.
PloS one volume16,
  pagese0253940 (year2021).

Beletal22
authorBelzile, L. R., authorDavison, A. C.,
  authorGampe, J., authorRootzén, H. &
  authorZholud, D.
journaltitleIs there a cap on longevity? A
  statistical review.
Annual Review of Statistics and Its
  Application volume9, pages21–45
  (year2022).

RooZho17
authorRootzén, H. & authorZholud, D.
journaltitleHuman life is unlimited–but
  short.
Extremes volume20,
  pages713–728 (year2017).

Baretal18
authorBarbi, E., authorLagona, F.,
  authorMarsili, M., authorVaupel, J. W. &
  authorWachter, K. W.
journaltitleThe plateau of human mortality:
  Demography of longevity pioneers.
Science volume360,
  pages1459–1461 (year2018).

Modetal17
authorModig, K., authorAndersson, T.,
  authorVaupel, J., authorRau, R. &
  authorAhlbom, A.
journaltitleHow long do centenarians survive?
  Life expectancy and maximum lifespan.
Journal of internal medicine
  volume282, pages156–163
  (year2017).

WilRob03
authorWilmoth, J. R. & authorRobine, J.-M.
journaltitleThe world trend in maximum life
  span.
Population and Development Review
  volume29, pages239–257
  (year2003).

GAvGav11
authorGavrilov, L. & authorGavrilova, N.
journaltitleMortality measurement at advanced
  ages: A study of the social security administration death master file.
North American Actuarial Journal
  volume15, pages442–447
  (year2011).

New18
authorNewman, S. J.
journaltitleErrors as a primary cause of
  late-life mortality deceleration and plateaus.
PLoS biology volume16,
  pagese2006776 (year2018).

Cam22
authorCamarda, C. G.
journaltitleThe curse of the plateau. measuring
  confidence in human mortality estimates at extreme ages.
Theoretical Population Biology
  volume144, pages24–36
  (year2022).

MisVau15
authorMissov, T. I. & authorVaupel, J. W.
journaltitleMortality implications of mortality
  plateaus.
siam REVIEW volume57,
  pages61–70 (year2015).

Gbaetal17
authorGbari, S., authorPoulain, M.,
  authorDal, L. & authorDenuit, M.
journaltitleExtreme value analysis of mortality
  at the oldest ages: a case study based on individual ages at death.
North American Actuarial Journal
  volume21, pages397–416
  (year2017).

HanSib16
authorHanayama, N. & authorSibuya, M.
journaltitleEstimating the upper limit of
  lifetime probability distribution, based on data of japanese centenarians.
Journals of Gerontology Series A: Biomedical
  Sciences and Medical Sciences volume71,
  pages1014–1021 (year2016).

Einetal19
authorEinmahl, J. J., authorEinmahl, J. H. &
  authorde Haan, L.
journaltitleLimits to human life span through
  extreme value theory.
Journal of the American Statistical
  Association volume114, pages1075–1080
  (year2019).

LiLiu20
authorLi, J. & authorLiu, J.
journaltitleA modified extreme value
  perspective on best-performance life expectancy.
Journal of Population Research
  volume37, pages345–375
  (year2020).

Mil20
authorMilholland, B.
journaltitleJeanne Calment, actuarial
  paradoxography and the limit to human lifespan.
Rejuvenation Research
  volume23, pages17–18
  (year2020).

tableurl
authorGroup, G. R.
titleWorld's oldest person titleholders since 1955.
howpublished<https://grg.org/Adams/C.HTM>
  (year2018).
note[Online; accessed 14-06-2023].

Bri86
authorBrillinger, D. R.
journaltitleA biometrics invited paper with
  discussion: the natural variability of vital rates and associated
  statistics.
Biometrics volume42,
  pages693–734 (year1986).

NM65
authorNelder, J. A. & authorMead, R.
journaltitleA Simplex Method for Function
  Minimization.
The Computer Journal
  volume7, pages308–313
  (year1965).

Kiss24
authorKiss, C.
titleModelling the age of the oldest person in the world.
howpublished<https://github.com/csabi0312/modelling-the-age-of-the-oldest-person-in-the-world>
  (year2024).

AbrSte84
authorAbramowitz, M. & authorStegun, I. A.
titlePocketbook of Mathematical Functions
  (publisherVerlag Harri Deutsch, addressThun-Frankfurt
  am Main, year1984).





§ ACKNOWLEDGEMENTS


We thank Katalin Kovács for some useful advice at an early stage of the project which led to this collaboration.
The work of Cs. Kiss and B. Vető was supported by the NKFI (National Research, Development and Innovation Office) grant FK142124.
B. Vető is also grateful for the support of the NKFI grant KKP144059 “Fractal geometry and applications” and for the Bolyai Research Scholarship of the Hungarian Academy of Sciences.
L. Németh was supported by MaRDI, funded by the Deutsche Forschungsgemeinschaft (DFG), project number 460135501, NFDI 29/1 “MaRDI – Mathematische Forschungsdateninitiative.



§ AUTHOR CONTRIBUTIONS


B. V. initiated the research. Cs. K. and B. V. derived the model, prepared the scripts and figures, and carried out the estimations. Cs. K., L. N. and B. V. analyzed the results and wrote the manuscript.



§ DATA AVAILABILITY


The titleholder data are freely available at https://grg.org/Adams/C.HTM



§ ADDITIONAL INFORMATION


The authors declare no competing interests.



§ TABLES








