Gompertz distribution
The Gompertz distribution is a continuous probability distribution defined on the positive real line, characterized by a probability density function of the form f(x; \eta, \theta) = \frac{\eta}{\theta} \exp\left( \frac{\eta x}{\theta} + 1 - \exp\left(\frac{\eta x}{\theta}\right) \right) for x > 0, where \eta > 0 is a shape parameter and \theta > 0 is a scale parameter, and featuring an exponentially increasing hazard rate that makes it suitable for modeling aging-related failure processes.[1] Named after the British actuary Benjamin Gompertz, the distribution derives from his 1825 proposal of a mathematical law describing human mortality as increasing geometrically with age, expressed through a force of mortality \mu(x) = B c^x for constants B > 0 and c > 1.[2] This law forms the basis of the distribution's survival function S(x) = \exp\left( -\int_0^x \mu(t) \, dt \right) = \exp\left( -a (c^x - 1) \right), where a = B / \ln c > 0, leading to the cumulative distribution function F(x) = 1 - \exp\left( -a (c^x - 1) \right).[1] Key properties include a hazard function h(x) = \frac{\eta}{\theta} \exp\left(\frac{\eta x}{\theta}\right) that starts low and rises exponentially, reflecting accelerating risk over time, and moments that involve the exponential integral function without closed-form expressions, such as the mean \mathbb{E}[X] = \frac{\theta}{\eta} \, e \, E_1(1), where E_1(z) = \int_z^\infty \frac{e^{-t}}{t} \, dt is the exponential integral.[3] The distribution is positively skewed, unimodal for certain parameter values, and related to extreme value distributions through transformations, enabling simulation via inverse CDF methods with uniform random variables.[1] The Gompertz distribution finds primary applications in actuarial science for pricing life insurance and annuities by modeling human mortality tables, in demography for analyzing population survival curves, and in reliability engineering for lifetimes of components subject to wear-out failures.[1] Extensions, such as the Gompertz-Makeham distribution incorporating a constant baseline hazard, further enhance its utility in epidemiological studies of age-specific death rates.[4]Definition
Parameterization
The Gompertz distribution is a continuous probability distribution for a nonnegative random variable X with support on [0, \infty). It is defined using two positive parameters: the scale parameter \eta > 0, and the shape parameter b > 0, which represents the rate of exponential increase in mortality. The parameter \eta appears in the cumulative hazard function as \eta (e^{b x} - 1), influencing the overall scale of the distribution, while b controls the acceleration of the hazard over time, capturing the aging process. This parameterization establishes the foundational hazard structure h(x) = b \eta e^{b x}, from which the probability density function and other distribution functions are derived. The initial hazard rate at age zero is h(0) = b \eta. Alternative parameterizations are employed in applications such as actuarial modeling of adult mortality, where a shifted form starts at an initial age x_0 > 0 and uses an adjusted initial rate c = b \eta e^{b x_0} to reflect the effective mortality at that age, while retaining the same exponential growth rate b. The standard unshifted form, however, remains the primary reference for theoretical developments and general survival analysis.Probability density function
The probability density function of the Gompertz distribution, parameterized by scale parameter \eta > 0 and shape parameter b > 0, is given by f(x; \eta, b) = b \eta \exp\left(\eta + b x - \eta e^{b x}\right), \quad x \geq 0. This form arises from the underlying Gompertz law of mortality, where the hazard rate h(x) = b \eta e^{b x} increases exponentially with age x. The survival function is then S(x) = \exp\left(-\int_0^x h(u) \, du\right) = \exp\left(-\eta (e^{b x} - 1)\right), and the PDF follows as f(x) = h(x) S(x), yielding the exponential structure above after simplification.[3] To verify normalization, note that \int_0^\infty f(x; \eta, b) \, dx = -\int_0^\infty \frac{d}{dx} S(x) \, dx = -[S(\infty) - S(0)] = 1 - 0 = 1, confirming it integrates to unity over [0, \infty).[3] The PDF shape is right-skewed, starting near zero at x = 0, rising to a single mode at x = \frac{1}{b} \ln\left(\frac{1}{\eta}\right) (provided \eta < 1; otherwise, the mode is at the boundary), and then decaying asymptotically to zero as x increases, reflecting accelerating failure rates typical in actuarial and reliability contexts.[5]Distribution functions
Cumulative distribution function
The cumulative distribution function (CDF) of the Gompertz distribution with parameters \eta > 0 and b > 0 is given by F(x; \eta, b) = 1 - \exp\left(-\eta (e^{b x} - 1)\right), \quad x \geq 0, and F(x; \eta, b) = 0 for x < 0. This formula arises as the integral of the probability density function from 0 to x, or equivalently as F(x) = 1 - S(x), where S(x) denotes the survival function. The CDF satisfies the boundary conditions F(0; \eta, b) = 0 and \lim_{x \to \infty} F(x; \eta, b) = 1. It is strictly increasing and continuous on [0, \infty), reflecting the positive support and non-negative density of the distribution. The quantile function, or inverse CDF, is the solution to F(q(p); \eta, b) = p for $0 < p < 1, yielding the explicit form q(p; \eta, b) = \frac{1}{b} \ln\left(1 + \frac{1}{\eta} \ln\left(\frac{1}{1 - p}\right)\right). This closed-form expression facilitates computation of quantiles, such as medians or percentiles, in applications like survival analysis.Survival function
The survival function of the Gompertz distribution, denoted S(x; \eta, b), represents the probability that a lifetime exceeds age x \geq 0, and is given by S(x; \eta, b) = \exp\left( -\eta \left( e^{b x} - 1 \right) \right), where \eta > 0 is a scale parameter and b > 0 is a shape parameter controlling the rate of increase in mortality.[6] This function is derived directly from the cumulative distribution function F(x; \eta, b) as S(x; \eta, b) = 1 - F(x; \eta, b), providing the complementary probability of survival beyond x after accounting for the accumulated risk up to that point.[6] Key properties include S(0; \eta, b) = 1, indicating certain survival at age zero, and the function is strictly decreasing to \lim_{x \to \infty} S(x; \eta, b) = 0, reflecting inevitable mortality over time.[6] Additionally, the natural logarithm of the survival function simplifies to \log S(x; \eta, b) = -\eta \left( e^{b x} - 1 \right), which highlights an exponential decay in survival probability, modulated by the age-dependent term e^{b x} that accelerates with increasing b.[6] In reliability and survival analysis, S(x; \eta, b) quantifies the likelihood of system or organism endurance past duration x, making it valuable for modeling aging processes with exponentially rising failure rates, such as human mortality in actuarial contexts.[6]Hazard rate function
The hazard rate function, also known as the force of mortality in actuarial science, for the Gompertz distribution is given by h(x; \eta, b) = b \eta e^{b x} for x \geq 0, where \eta > 0 is a scale parameter and b > 0 is a shape parameter controlling the rate of increase in mortality.[7] This explicit form is derived from the general definition of the hazard rate as the ratio of the probability density function to the survival function, h(x) = f(x)/S(x), or equivalently as h(x) = -\frac{d}{dx} \log S(x). Applying these to the Gompertz distribution's component functions results in the characteristic linear exponential expression.[7] Key properties of the hazard function include its strict monotonicity: it is increasing for all x \geq 0 because the derivative h'(x) = b^2 \eta e^{b x} > 0. The initial value is h(0) = b \eta, and \lim_{x \to \infty} h(x) = \infty, reflecting an unbounded escalation in failure risk over time that suits modeling accelerating mortality.[7] This functional form captures the Gompertz law of mortality, first articulated by Benjamin Gompertz in 1825, which describes how the intensity of mortality increases exponentially with age owing to a progressive, uniform weakening of the body's capacity to resist death across successive equal time periods.[8][9]Moments and generating functions
Raw moments
The raw moments of the Gompertz distribution, assuming the hazard function h(x) = a e^{b x} for x \geq 0 with parameters a > 0 and b > 0, are given by the general formula E[X^r] = \frac{r!}{b^r} \, e^{a/b} \, E_{r,1}\left( \frac{a}{b} \right), where E_{n,s}(z) = \int_1^\infty e^{-z t} t^{-n} \, dt is the generalized exponential integral function. The first raw moment, or mean, is thus E[X] = \frac{1}{b} \, e^{a/b} \, E_1\left( \frac{a}{b} \right), with E_1(z) = \int_z^\infty \frac{e^{-t}}{t} \, dt denoting the standard exponential integral; equivalently, E[X] = -\frac{1}{b} \, e^{a/b} \, \Ei\left( -\frac{a}{b} \right), where \Ei is the exponential integral function. This expression quantifies the central tendency, reflecting the distribution's heavy-tailed nature influenced by the aging rate b and initial hazard level a. The second raw moment is E[X^2] = \frac{2}{b^2} \, e^{a/b} \, E_2\left( \frac{a}{b} \right), where E_2(z) follows from the generalized form. The variance is obtained as \Var(X) = E[X^2] - (E[X])^2, yielding the explicit relation \Var(X) = \frac{2}{b^2} \, e^{a/b} \left[ -\frac{a}{b} \ {}_3F_3\left(1,1,1;2,2,2;\frac{a}{b}\right) + \frac{1}{2} \left( \frac{\pi^2}{6} + \left( \gamma + \ln \frac{a}{b} \right)^2 \right) \right] - (E[X])^2, with \gamma \approx 0.5772156649 the Euler-Mascheroni constant and {}_3F_3 the generalized hypergeometric function; this measures the spread, approaching \pi^2/(6 b^2) for small a/b. Higher-order raw moments follow the general formula involving E_{r,1}, but practical computation is challenging due to the need for numerical evaluation of the special functions, particularly for small a/b where singularities and logarithmic terms complicate stability.Moment-generating function
The moment-generating function (MGF) of the Gompertz distribution, denoted M(t) = \mathbb{E}[e^{tX}], where X follows the distribution with parameters b > 0 and \eta > 0, is derived by integrating the probability density function: M(t) = \int_0^\infty e^{tx} \, b \eta e^{bx} \exp\left(-\eta (e^{bx} - 1)\right) \, dx. This integral can be evaluated using the substitution u = e^{bx}, which yields dx = du / (b u) and transforms the expression into a form involving the upper incomplete gamma function \Gamma(s, z) = \int_z^\infty v^{s-1} e^{-v} \, dv: M(t) = e^\eta \eta^{-t/b} \Gamma\left(1 + \frac{t}{b}, \eta\right), valid for all real t since all moments exist due to the rapid decay of the distribution's tail.[3] The derivation proceeds by simplifying the exponent in the integrand to (t + b)x - \eta e^{bx}, applying the substitution to obtain \eta e^\eta \int_1^\infty u^{t/b} e^{-\eta u} \, du, and then rescaling the variable v = \eta u to recognize the incomplete gamma form after adjusting the power and limits. This relation to the incomplete gamma function also connects the MGF to the Laplace transform of the distribution, which is M(-t).[3] The MGF facilitates the computation of moments and cumulants through differentiation: the k-th raw moment is M^{(k)}(0). However, explicit evaluation of M(t) or its derivatives often requires numerical methods, as the incomplete gamma function lacks a simpler closed-form expression in elementary functions for arbitrary parameters. The raw moments, obtainable via this differentiation, provide key descriptive statistics but are not expressible in elementary terms for general k.Parameter estimation
Maximum likelihood estimation
Maximum likelihood estimation (MLE) provides optimal estimators for the parameters η and b of the Gompertz distribution under standard regularity conditions, utilizing the full information from the observed data.[3] This approach is particularly valuable in applications involving lifetime data, where the distribution often arises in modeling increasing failure rates.[10] For a random sample of n uncensored observations x_1, \dots, x_n from the Gompertz distribution, the likelihood function is given by L(\eta, b \mid \mathbf{x}) = \prod_{i=1}^n f(x_i \mid \eta, b), where f(x \mid \eta, b) is the probability density function.[3] In the presence of right-censoring, the likelihood incorporates the survival function S(x \mid \eta, b) for censored observations, yielding L(\eta, b \mid \mathbf{x}, \boldsymbol{\delta}) = \prod_{i=1}^n \left[ f(x_i \mid \eta, b) \right]^{\delta_i} \left[ S(x_i \mid \eta, b) \right]^{1 - \delta_i}, with \delta_i = 1 for uncensored and \delta_i = 0 for censored cases.[3] The log-likelihood for uncensored data simplifies to \begin{align*} l(\eta, b \mid \mathbf{x}) &= \sum_{i=1}^n \left[ \log(b \eta) + \eta + b x_i - \eta e^{b x_i} \right] \ &= n \log(b \eta) + n \eta + b \sum_{i=1}^n x_i - \eta \sum_{i=1}^n e^{b x_i}. \end{align*} Maximizing this function requires solving the score equations \partial l / \partial \eta = 0 and \partial l / \partial b = 0, which generally lack closed-form solutions due to their nonlinearity and must be addressed via numerical optimization methods such as Newton-Raphson iteration.[3] For censored data, the log-likelihood takes a similar form but adjusts the summation over uncensored terms for the density contributions.[10] Under standard assumptions, the MLEs \hat{\eta} and \hat{b} are consistent and asymptotically efficient, with asymptotic normality \sqrt{n} (\hat{\theta} - \theta) \to \mathcal{N}(0, \mathcal{I}(\theta)^{-1}), where \theta = (\eta, b) and \mathcal{I}(\theta) is the Fisher information matrix.[10] Standard errors for inference can be obtained from the inverse of the observed Hessian matrix evaluated at the MLE.[3] The nonlinear nature of the optimization can lead to convergence difficulties, especially with small sample sizes or highly skewed data, potentially requiring careful initial parameter guesses or alternative algorithms like BFGS.[10] Implementations are available in statistical software, such as themlgompertz function in the R package univariateML for uncensored data using Newton-Raphson.[11] For survival analysis with censoring, R's flexsurv package supports Gompertz models via parametric fitting.
Method of moments estimation
The method of moments (MoM) estimation for the two-parameter Gompertz distribution equates the sample mean \bar{x} and sample variance s^2 to the theoretical mean E[X] and variance \mathrm{Var}(X), respectively, yielding a system of nonlinear equations solved numerically for the parameters \eta > 0 and b > 0.[3] The theoretical mean is given by E[X] = \frac{e^{\eta} E_1(\eta)}{b}, where E_1(\cdot) denotes the exponential integral function E_1(z) = \int_z^\infty \frac{e^{-t}}{t} \, dt. The variance is \mathrm{Var}(X) = \frac{2}{b^2} e^{\eta} \left[ -\eta \, {}_3F_3\left(1,1,1; 2,2,2; \eta \right) + \frac{1}{2} \left( \frac{\pi^2}{6} + (\gamma + \ln \eta)^2 \right) \right] - (E[X])^2, with \gamma \approx 0.57721 the Euler-Mascheroni constant and {}_3F_3(\cdot) the generalized hypergeometric function; alternatively, it can be computed from the second raw moment minus the squared mean, where raw moments are derived via the Laplace transform or integration.[3] To solve, an initial estimate of b is obtained from the variance-mean relation, often using the approximation \mathrm{Var}(X) \approx \pi^2 / (6 b^2) valid for small \eta, so b \approx \pi / \sqrt{6 s^2}, though the full relation requires iteration to account for \eta-dependence. With b fixed, \eta is then found by numerically inverting e^{\eta} E_1(\eta) = b \bar{x} via methods like Newton-Raphson, as no closed-form solution exists.[3] This approach leverages the raw moments of the distribution for straightforward computation from sample summaries. MoM offers simplicity and serves as a quick initial estimate for more sophisticated methods, but it is generally less statistically efficient than maximum likelihood estimation and can be sensitive to deviations from the assumed distributional form, particularly in the tails.[12]Characteristic properties
Parameter effects on shape
In the general parameterization of the Gompertz distribution with hazard function h(x) = \eta \exp(b x) for \eta > 0 (initial hazard rate) and b > 0 (growth rate, independent of \eta), larger values of \eta elevate the baseline hazard across the support, resulting in higher initial probability density and a leftward shift in the mass of the distribution toward lower values of x, as failures become more likely at earlier stages. This effect is evident in the probability density function (PDF), where increased \eta amplifies the density near x = 0 while compressing the tail, thereby raising the early hazard rate and reducing the expected lifetime.[13] In contrast, larger values of b accelerate the exponential growth of the hazard function, leading to a slower initial accumulation of risk followed by a more rapid escalation later on. This produces a more pronounced right skew in the PDF and cumulative distribution function (CDF), with the bulk of the probability mass concentrating at higher ages or times, reflecting delayed but intensified failure risks. The survival function decays more gradually at first but drops sharply thereafter, emphasizing the distribution's utility in modeling accelerating aging processes.[3] Visualizations of parameter families illustrate these dynamics clearly. For fixed b, plots of the PDF with varying \eta show densities starting higher and tapering more quickly as \eta increases, often transitioning from unimodal to monotonically decreasing shapes when \eta exceeds b. Similarly, for fixed \eta, increasing b yields families of PDFs and hazard rates that exhibit delayed peaking, with the mode shifting rightward and the hazard curve steepening exponentially; the CDF rises more slowly initially before accelerating. The mode, when it exists (for b > \eta), is located at x = \frac{1}{b} \ln\left(\frac{b}{\eta}\right), though this position sensitively depends on the relative magnitudes of the parameters.[7] Small perturbations in b are particularly influential in fitting real-world data, as they capture subtle accelerations in aging or failure rates, such as in human mortality tables where even minor increases in b align the model with observed rises in late-life hazards without overpredicting early events. This parameter sensitivity underscores b's role in tailoring the distribution to empirical trends in demography and reliability.[3]Information measures
The Kullback-Leibler (KL) divergence is a measure of the difference between two probability distributions, defined as D_{\text{KL}}(P \parallel Q) = \int f_P(x) \log \left( \frac{f_P(x)}{f_Q(x)} \right) dx, where f_P and f_Q are the probability density functions of distributions P and Q. For the Gompertz distribution, a closed-form expression exists for the KL divergence between two Gompertz distributions with parameters (b_1, q_1) and (b_2, q_2), where the pdf is given by f(x \mid b, q) = e^q b q e^{b x} e^{-q e^{b x}} for x \geq 0 and b, q > 0. This expression is D_{\text{KL}}(F_1 \parallel F_2) = \ln\left( \frac{e^{q_1} b_1 q_1}{e^{q_2} b_2 q_2} \right) + e^{q_1} \left[ \left( \frac{b_2}{b_1} - 1 \right) \text{Ei}(-q_1) + \frac{q_2}{q_1} q_1^{-b_2/b_1} \Gamma\left( \frac{b_2}{b_1} + 1, q_1 \right) \right] - (q_1 + 1), involving the exponential integral \text{Ei}(\cdot) and the upper incomplete gamma function \Gamma(s, x).[14] Computing the KL divergence for a Gompertz distribution against alternatives like the Weibull distribution typically requires numerical integration due to the lack of closed forms, as the integral involves products of exponentials that do not simplify analytically. Approximations or quadrature methods, such as Gauss-Laguerre integration, are commonly employed for such cases. Exact expressions for Gompertz-Weibull comparisons remain unavailable in closed form, though general derivations for univariate continuous distributions provide frameworks using moment-generating functions for numerical evaluation.[15] In applications, the KL divergence aids model selection by quantifying how well a Gompertz model approximates empirical data compared to alternatives, and supports goodness-of-fit testing in lifetime analysis, such as mortality modeling where distinguishing Gompertz from Weibull fits is crucial.[14][15] Other information measures for the Gompertz distribution include the Shannon entropy, defined as H = -\int f(x) \log f(x) \, dx, which quantifies uncertainty in the distribution. No closed-form expression exists for the Shannon entropy of the standard Gompertz distribution, necessitating numerical estimation via methods like Monte Carlo integration or series expansions. Comparative studies of multiple entropy measures (e.g., Rényi, Tsallis) for Gompertz distributions also rely on numerical computation to assess properties like truncation effects in reliability contexts.[16]Related distributions
Gompertz–Makeham distribution
The Gompertz–Makeham distribution extends the Gompertz distribution by incorporating an additional constant term in the hazard function to account for age-independent mortality risks, such as accidents. The hazard function is given byh(t) = a + b e^{c t},
where t \geq 0 represents age, a > 0 is the constant baseline hazard, b > 0 scales the initial level of the age-dependent component, and c > 0 governs the exponential rate of increase with age.[17] This formulation was originally proposed by William Matthew Makeham to improve the modeling of human mortality patterns observed in actuarial data.[18] The probability density function (PDF) and cumulative distribution function (CDF) derive from the hazard via the survival function S(t) = \exp\left( -\int_0^t h(u) \, du \right), yielding
f(t) = (a + b e^{c t}) \exp\left( -a t - \frac{b}{c} (e^{c t} - 1) \right)
for the PDF and
F(t) = 1 - \exp\left( -a t - \frac{b}{c} (e^{c t} - 1) \right)
for the CDF.[17] These closed-form expressions facilitate analytical computations in survival analysis without requiring numerical integration. When the baseline parameter a = 0, the distribution reduces to the standard Gompertz distribution with hazard b e^{c t} and parameters b and c.[17] The inclusion of the positive constant a captures mortality causes unrelated to aging, such as external hazards, which enhances the model's empirical fit to human lifetime data across various populations by addressing deviations from pure exponential growth in observed death rates at younger ages.[17]