A novel three-parameter Gompertz-Lomax Distribution with a physics-informed survival network: statistical theory, hazard shape characterization, and real-world applications

1 Introduction

Precise statistical modeling of lifetime data is fundamental to reliability engineering, clinical survival analysis, and risk management across virtually all applied sciences. A persistent challenge is that real-world failure processes exhibit qualitatively different hazard rate behaviors that no single two-parameter distribution accommodates simultaneously: decreasing hazard during early-life burn-in, increasing hazard during aging and wear-out, and the bathtub shape that combines both (Murthy et al., 2004).

The Gompertz distribution (Gompertz, 1825) captures monotone increasing hazard, but cannot model decreasing or bathtub shapes. The Lomax distribution (Lomax, 1954) is restricted to decreasing hazard. The Weibull distribution (Weibull, 1951) covers increasing and decreasing, but not the bathtub shape (Mudholkar et al., 1996). Three-parameter extensions, such as the Modified Weibull (Lai et al., 2003) and the Flexible Weibull (Bebbington et al., 2007) have addressed the bathtub requirement, but typically at the cost of closed-form tractability.

A complementary strand of research has produced analytically elegant one- and two-parameter distributions through mixture and transformation constructions. Notable contributions include the Shanker distribution (Shanker, 2015), the Sujatha distribution (Shanker, 2016), the Rani distribution (Shanker, 2017), the Ram Awadh distribution (Shukla, 2018), and several others applied to medical, reliability, and engineering data. Extensions using exponentiation and discretization have been developed for count data and uncertain observations in these domains. However, the mixture-based one-parameter constructions in this family are structurally limited to decreasing or unimodal hazard profiles and cannot represent the bathtub shape, critical for many industrial and biomedical applications.

Three gaps persist simultaneously in the literature. First, the natural first-principles approach to a bathtub distribution, namely, directly summing a Gompertz increasing hazard component with a Lomax decreasing hazard component, has not been explored as a standalone three-parameter family. Second, no deep learning survival model simultaneously outputs the parameters of an analytically tractable three-parameter distribution and enforces the hazard shape conditions that follow from distributional theory (Wang et al., 2024; Klambauer et al., 2020). Third, goodness-of-fit evaluation in distribution articles seldom reports all standard information criteria alongside Q-Q plots and frequency diagrams in a single paper. No existing three-parameter model provides a closed-form survival function that analytically separates IFR and bathtub regimes through a single testable algebraic condition on its parameters. Recent extensions of the Lomax family have explored flexible hazard shapes (Ul Hassan and Dar, 2020; Al-Babtain et al., 2021), and Bayesian estimation under progressive censoring for Gompertz-type models has been studied by Almetwally et al. (2022); however, none of these works combines the additive composite hazard structure with a physics-constrained deep learning component.

1.1 Novelty relative to existing three-parameter lifetime models

Several three-parameter lifetime models exist in the literature, and the novelty of the Gompertz-Lomax Distribution (GLD) must be assessed against this body of work. The Modified Weibull (Lai et al., 2003) and Flexible Weibull (Bebbington et al., 2007) achieve bathtub hazard but do not admit a closed-form survival function free of exponential integrals. The Exponentiated Gompertz (El-Damcese et al., 2015) produces an increasing or unimodal hazard only. The Kumaraswamy distribution (Jones, 2009) has a flexible density but is bounded and lacks a bathtub hazard. The generalized Lomax families (Ul Hassan and Dar, 2020; Al-Babtain et al., 2021) produce bathtub hazard but not through an additive composite construction. The key distinguishing features of the GLD are: (i) its hazard is the additive sum of an IFR component and a DFR component, providing a direct first-principles motivation; (ii) the survival function is fully elementary and closed-form (no special functions required for S, f, F, or h); and (iii) the bathtub condition ab≷c² is a single, testable algebraic inequality on the estimated parameters, a property not shared by other three-parameter bathtub models. These three features together justify the GLD as a distinct contribution rather than a marginal extension of an existing family.

The primary objective of this work is to construct and fully characterize such a distribution, derive its complete analytical properties, develop a physics-constrained deep survival network that outputs its parameters, and validates both contributions empirically across multiple application domains.

The GLD has hazard

h(x;a,b,c)=aebx+c1+cx, x>0,a,b,c>0.(1)

The Gompertz term ae^*bx* is monotone increasing, while the Lomax term c/(1+cx) is monotone decreasing. Their sum can therefore produce increasing, decreasing, or bathtub-shaped hazard depending on whether the Gompertz component dominates in the long run or the Lomax component dominates initially. The resulting survival function is fully closed-form (Section 2), and the hazard shape theorem (Theorem 3.2) provides an exact, testable condition separating the three regimes.

The remainder of this article is organized as follows. Section 2 defines the GLD and establishes its distributional construction and special cases. Section 3 derives the complete statistical and reliability properties. Section 4 presents the four estimation procedures. Section 5 reports the Monte Carlo simulation study, including a comparative evaluation of all four estimators. Section 6 introduces the Physics-Informed Survival Network (PISN), which embeds the distributional hazard shape constraint as a differentiable training penalty in a deep neural network. Section 7 applies both the GLD and PISN to three benchmark datasets. Section 8 concludes with a summary of contributions, limitations, and future directions.

2 The Gompertz-Lomax Distribution

2.1 Construction

Integrating h(x) from Equation 1 yields the cumulative hazard H(x)=∫0xh(t)dt=(a/b)(ebx-1)+ln(1+cx), from which the survival, density, and distribution functions follow via S(x) = e^−*H*(*x*) and f(x) = h(x)S(x). It is important to note that while S(x), f(x), F(x), and h(x) are all analytically closed-form (explicit elementary functions of x and the parameters), some derived quantities such as moments, the quantile function, and the mean residual life do not admit simple closed forms and are evaluated numerically; this distinction is made explicit throughout Section 3.

Definition 0.0.1. X~GLD(a, b, c) if its survival, density, and distribution functions are

S(x;a,b,c)=exp(-ab(ebx-1))(1+cx)-1,(2)

f(x;a,b,c)=(aebx+c1+cx)S(x;a,b,c),(3)

F(x;a,b,c)=1-S(x;a,b,c),(4)

for x>0, a, b, c>0.

Verification that ∫0∞fdx=1 follows from [-S]0∞=1, since S(0) = 1 and S(∞) = 0. The distribution function (Equation 4) follows directly from the survival function.

Proposition 2.2 (Special cases).

• (i) GLD(a, b, 0)≡Gompertz(a, b).
• (ii) As a → 0, S(x) → (1+cx)⁻¹, so the GLD converges to a Lomax distribution with shape α = 1 and scale λ = 1/c, i.e., Lomax(α = 1, λ = 1/c).
• (iii) As b → 0 with a/b → μ, it converges to a modified exponential-Lomax model.
• (iv) As a → λ, b → 0, c → 0, the distribution converges to Exponential(λ).

Figures 1, 2 display the PDF and CDF for representative parameter combinations.

3 Statistical and reliability properties

3.1 Moments

Theorem 3.1 (Raw moments). For integer r≥1, μr′=E[Xr] is given by the convergent series (Equation 5).

μr′=ea/b∑k=0∞(-a/b)kk!∫0∞xr(aebx+c1+cx)e-kebx(1+cx)-1dx,(5)

where each integral reduces to a combination of the exponential integral E₁(·) via the substitution u = e^*bx*.

Proof. Write μr′=∫0∞xrf(x)dx=∫0∞xrh(x)S(x)dx. Substitute the closed-form survival function S(x) = exp(−(a/b)(e^*bx*−1))(1+cx)⁻¹ and expand the Gompertz exponential factor as exp(-(a/b)(ebx-1))=ea/b∑k=0∞(-a/b)kk!e-kebx, which converges absolutely for all a/b>0 and all x≥0. Interchange of summation and integration is justified by the dominated convergence theorem (the partial sums are bounded by an integrable function for a/b ≤ 10). Each resulting integral ∫0∞xre-kebx(1+cx)-1dx is then reduced via the substitution u = ke^*bx* to a finite linear combination of exponential integrals E₁(·) and incomplete gamma functions, all of which are computable in closed form. In practice, Gauss-Laguerre quadrature is used for numerical stability.

In practice, μr′ is computed numerically via Gauss-Laguerre quadrature; the series converges rapidly for a/b ≤ 10. Central moments μr=E[(X-μ1′)r] follow by the standard recurrence. The coefficient of variation, skewness β1=μ3/μ23/2, kurtosis β2=μ4/μ22, and dispersion index γ=μ2/μ1′ are plotted in Figure 3.

Following Khan et al. (2021), a CDF-based skewness measure is defined as

SK=2F(μ′)-1,(6)

where μ′ = E[X] is the distribution mean and F(·) is the cumulative distribution function. For a perfectly symmetric distribution F(μ′) = 0.5 so S_*K* = 0; positive (negative) values indicate right (left) skewness. Unlike the conventional moment-based coefficient β1=μ3/μ23/2, the measure S_*K* depends only on the CDF evaluated at a single point and is therefore robust to heavy tails and free from higher-order moment existence conditions (Khan et al., 2021).

Substituting the GLD cumulative distribution function F(x; a, b, c) = 1 − exp((a/b)(1−e^*bx*))(1+cx)⁻¹ into Equation 6 yields the expression

SK(a,b,c)=2[1-e(a/b)(1-ebμ′)(1+cμ′)-1]-1,(7)

where μ′ is the GLD mean evaluated numerically via Gauss–Laguerre quadrature as described above. Expression Equation 7 is closed-form in μ′ and is straightforward to evaluate for any given parameter triplet (a, b, c).

The behavior of S_*K* across the GLD parameter space is illustrated in Figure 4. Panels (a) and (b) show S_*K* as a function of c for varying a and b respectively (remaining parameter fixed). Panel (c) shows S_*K* as a function of a for varying c. Panel (d) presents a joint heatmap of S_*K* over the (a, c)-plane with b = 0.5 fixed, together with the hazard-shape boundary c=ab (dashed) from Theorem 0.0.4. Several features are noteworthy. First, S_*K*>0 throughout most of the parameter space, confirming that the GLD is predominantly right-skewed, consistent with lifetime applications. Second, in the bathtub regime (ab<c²), the Lomax-dominated phase reduces S_*K* and can produce near-zero or mildly negative skewness. Third, increasing a monotonically increases S_*K* for fixed b and c, while increasing c generally decreases S_*K*, reflecting the competing contributions of the two hazard components. The moment-based skewness β1, kurtosis β₂, coefficient of variation, and dispersion index are plotted in Figure 3.

3.2 Incomplete moments and mean deviations

The s-th incomplete moment ϕs(t0)=∫0t0xsf(x)dx is computed using the same series as Theorem 0.0.3 with upper integration limit t₀. The mean deviations from the mean and median are δ1=2μ1′F(μ1′)-2ϕ1(μ1′) and δ2=μ1′-2ϕ1(M) respectively, where M denotes the median (Shaked and Shanthikumar, 2007).

3.3 Quantile function and median

The quantile function Q(p) = F⁻¹(p) satisfies the implicit equation

ab(ebQ-1)+ln(1+cQ)=-ln(1-p),(8)

obtained by inverting F(x) = 1−S(x). This equation does not have an elementary closed-form solution for Q(p), and is solved numerically via Brent's method (Brent, 1973) using the Gompertz quantile as the initial value. The median is M = Q(0.5).

3.4 Hazard rate and shape analysis

The hazard rate (Equation 1) is displayed in Figure 5 for IFR, DFR, and bathtub configurations.

Theorem 3.2 (Hazard Shape Characterization). Let h′(x) = abe^*bx*−c²/(1+cx)².

• *If ab>c²: h′(x)>0 for all x>0, so h is strictly IFR.*
• *If ab<c²: h′(0) < 0 and h′(x) → +∞ as x → ∞, so a unique turning point x^* exists satisfying abe^*bx** = c²/(1+c^{x^*)2}, giving a bathtub shape (UBT).*
• *If a → 0: h′(x) < 0 everywhere, so h is DFR (Lomax limit).*

Proof. Differentiate h(x) = ae^*bx*+c/(1+cx) to obtain h′(x) = abe^*bx*−c²/(1+cx)². At x = 0: h′(0) = ab−c². The first term abe^*bx* is strictly increasing and unbounded, while the second term c²/(1+cx)² is strictly decreasing to zero. Hence h′ has at most one zero. If ab>c², then h′(0)>0 and since the first term dominates for all x>0, we have h′>0 everywhere, giving IFR. If ab<c², then h′(0) < 0, but h′(x) → +∞, so by continuity there exists a unique x^*>0 where h′(x^*) = 0, confirming the bathtub shape.

Theorem 3.3 (Stochastic ordering). Let X_*i*~GLD(a_*i*, b, c), i = 1, 2. If a₁≥a₂ then X₁ ≤_lrX₂, implying hazard rate, mean residual life, and usual stochastic orders.

Proof. The likelihood ratio is

f1(x)f2(x)=a1ebx+c/(1+cx)a2ebx+c/(1+cx)·exp(-a1-a2b(ebx-1)).

Denote the first factor by R(x) and the second by E(x). Since a₁>a₂, we have E(x) strictly decreasing in x (its exponent is negative and monotone decreasing). To show R(x) is also strictly decreasing, differentiate:

R′(x)=(a1-a2)bebx·c/(1+cx)2-0[a2ebx+c/(1+cx)]2·(-1)=-(a1-a2)bebxc(1+cx)-2[a2ebx+c/(1+cx)]2.

Because a₁>a₂, b>0, and c>0, every factor in the numerator is strictly positive, so R′(x) < 0 for all x>0. Therefore f₁(x)/f₂(x) is the product of two strictly decreasing positive functions, hence strictly decreasing, establishing X₁ ≤_lrX₂. The chain lr ⇒ hr ⇒ mrl ⇒ st follows from standard theory (Shaked and Shanthikumar, 2007).

Remark 3.4. The condition ab≷c² in Theorem 0.0.4 provides a direct, testable diagnostic from the MLE: the fitted parameter triplet (â,b^,ĉ) immediately identifies the hazard regime of the population. This condition is also the physics constraint encoded in the PISN loss (Section 6).

3.5 Mean residual life

The mean residual life at age x₀ is

m(x0)=1S(x0)∫x0∞S(t)dt,(9)

so m(x₀) (Equation 9) is evaluated numerically via Gauss-Laguerre quadrature. Where S(t) is the closed-form survival function (Equation 2). Although S(t) is explicit, the integral ∫x0∞S(t)dt does not reduce to elementary functions, so m(x₀) is evaluated numerically via Gauss-Laguerre quadrature and displayed in Figure 6. For the IFR regime (ab>c²), m(x₀) is strictly decreasing. For the bathtub regime (ab<c²), m(x₀) initially increases before decreasing, consistent with early-life reliability improvement followed by long-run aging. The identity h(x) = [1+m′(x)]/m(x) provides an independent verification of Theorem 0.0.4.

3.6 Stress-strength reliability

Let X~GLD(a₁, b, c₁) denote component strength and Y~GLD(a₂, b, c₂) denote operational stress, sharing shape parameter b. The stress-strength reliability is

R=P(Y<X)=1-e(a1+a2)/b∑j,k=0∞(-a1/b)j(-a2/b)kj!k!Kj,k,(10)

The stress-strength reliability (Equation 10) is evaluated via the substitution u = e^*bx*. Where Kj,k=∫0∞[a1ebx+c1/(1+c1x)]e-(j+k)ebx(1+c2x)-1dx is evaluated by the substitution u = e^*bx* and reduces to a combination of incomplete gamma functions (Shaked and Shanthikumar, 2007). The MLE of R follows by invariance: R^=R(â1,b^,ĉ1,â2,ĉ2).

3.7 Order statistics

The PDF of the k-th order statistic from an i.i.d. sample of size n is

fX(k)(x)=n!(k-1)!(n-k)![F(x)]k-1[S(x)]n-kf(x).(11)

Expanding the order statistic density (Equation 11) using the binomial theorem [F(x)]^*k*−1 = (1−S(x))^*k*−1 via the binomial theorem and substituting Equations 2, 3 yields a closed expression involving GLD survival functions evaluated at stretched Gompertz parameters a(n−k+j+1).

3.8 Rényi entropy

For δ>0, δ≠1, the Rényi entropy is TR(δ)=(1-δ)-1ln(∫0∞f(x)δdx). Using f(x) = h(x)S(x) and S(x)^δ = exp(−δh(x)), where h(x) = (a/b)(e^*bx*−1)+ln(1+cx), we obtain

∫0∞f(x)δdx=eδa/b∑k=0∞(-δa/b)kk!∫0∞[h(x)]δe-kebx(1+cx)-δdx,(12)

where the expansion e-δ(a/b)(ebx-1)=eδa/b∑k=0∞(-δa/b)ke-kebx/k! is used, and each inner integral is evaluated via the substitution u = e^*bx*−1 combined with the binomial expansion of [h(x)]^δ = (a^{e^*bx*+c/(1+cx))δ}. Substituting Equation 12 into the entropy definition yields.

TR(δ)=11−δln​{eδa/b∑k=0∞(−δa/b)kk!∫0∞[h(x)]δe−k(ebx−1)(1+cx)−δdx}.(13)

As δ → 1 Equation 13 converges to the Shannon entropy entropy h=-∫0∞fln fdx by L'hôpital's rule.

3.9 Bonferroni and Lorenz curves

The Lorenz and Bonferroni curves are

L(p)=ϕ1(Q(p))μ1′, B(p)=L(p)p,(14)

where ϕ₁ is the first incomplete moment and Q(p) is the quantile function (Equation 8). The Gini coefficient G=1-2∫01L(p)dp provides a scalar summary of lifetime inequality. The Lorenz and Bonferroni curves (Equation 14) are plotted in Figure 7. Figure 7 shows L(p) and B(p) for three parameter configurations.

4 Parameter estimation

4.1 Maximum likelihood estimation

For an i.i.d. sample (x₁, …, x_*n*), the log-likelihood (Equation 15) is maximized by solving Equations 16–18.

ℓ(a,b,c)=∑i=1nln(aebxi+c1+cxi)-ab∑i=1n(ebxi-1)

-∑i=1nln(1+cxi).(15)

The score equations ∂ℓ/∂a = 0, ∂ℓ/∂b = 0, ∂ℓ/∂c = 0 are The MLE θ^=(â,b^,ĉ)⊤ is obtained by solving Equations 16–18 via L-BFGS-B (Byrd et al., 1995) with analytical gradients, initialized by differential evolution (Storn and Price, 1997) for global exploration. Standard errors are SE(θ^j)=[I-1(θ^)]jj1/2, where I is the observed Fisher information matrix. Standard errors are reported alongside parameter estimates in Tables 1–3.

∂ℓ∂a=∑i=1nebxiaebxi+c/(1+cxi)-1b∑i=1n(ebxi-1)=0,(16)

∂ℓ∂b=∑i=1naxiebxiaebxi+c/(1+cxi)+ab2∑i=1n(ebxi−1)−ab∑i=1nxiebxi=0,(17)

∂ℓ∂c=∑i=1n(1+cxi)-2aebxi+c/(1+cxi)-∑i=1nxi1+cxi=0.(18)

Table 1. Goodness-of-fit statistics for DS1 (leukemia survival, n = 33).

MLE estimates with standard errors (SE) in parentheses. Best values per criterion in bold. ^‡See footnote for the Modified Weibull (MW) entry.

The MW fit yielded an anomalously large K-S statistic on DS1 and DS2, attributable to the optimizer converging to a local minimum of the MW likelihood for these datasets. Multiple restarts with differential evolution and grid initialization was attempted; the result was reported reflects the best likelihood found across all starting values. The MW is a capable model on other lifetime datasets; the poor fit here is dataset-specific.

Table 2. Goodness-of-fit statistics for DS2 (COVID-19 survival, n = 50).

MLE estimates with standard errors (SE) in parentheses. Best values per criterion in bold.

See footnote to Table 1 regarding the MW anomalous fit.

Table 3. Goodness-of-fit statistics for DS3 [electronic component failures, n = 50; (Aarset, 1987)].

MLE estimates with standard errors (SE) in parentheses. Best values per criterion in bold.

4.2 Method of moments, LSE, and WLSE

The method of moments equates μr′(a,b,c)=n-1∑ixir for r = 1, 2, 3 and solves via Newton-Raphson. The least squares estimator (LSE) minimizes ∑i[F(x(i))-i/(n+1)]2 and the weighted LSE (WLSE) uses weights wi=(n+1)2(n+2)/[i(n-i+1)], both solved by L-BFGS-B (Byrd et al., 1995). The finite-sample performance of all four estimators is compared in the simulation study (Section 5).

5 Simulation study

The finite-sample performance of all four estimators (MLE, MOM, LSE, and WLSE) was assessed through N = 10, 000 Monte Carlo replications at n∈{25, 50, 100, 200, 300} for four parameter configurations: Set I (0.05, 0.30, 1.00), Set II (0.10, 0.20, 0.50), Set III (0.50, 0.50, 0.30), Set IV (0.10, 0.50, 0.50). Random samples were generated via the inverse transform method with Brent's quantile solver. Performance metrics are average bias (AB) and mean squared error (MSE):

AB(θ^)=1N∑i=1N(θ^i-θ), MSE(θ^)=1N∑i=1N(θ^i-θ)2.(19)

AB and MSE (Equation 19) are reported in Table 4, with AB and MSE decreasing monotonically with n for all four parameter sets, confirming MLE consistency. Table 5 provides a comparative summary of AB and MSE for all four estimators at representative sample sizes. Lower MSE is preferred; MLE achieves the lowest MSE at moderate to large n, while MOM and WLSE offer competitive performance at small n where parameter interactions can cause MLE to be sensitive to initialization. Figure 8 presents MSE heatmaps for visual comparison. Figure 9 shows MSE line profiles for Sets I and II.

Table 4. Average bias (AB) and mean squared error (MSE) of the GLD MLE under four parameter configurations and five sample sizes (N = 10, 000 replications).

Table 5. Comparative MSE of MLE, MOM, LSE, and WLSE for parameter ĉ under Set I and Set III at selected sample sizes (N = 10, 000 replications).

6 Physics-informed survival network

6.1 Architecture and forward pass

The forward pass (Equations 20, 21) maps covariates to GLD parameters. Let xi∈ℝp be covariates for subject i, and (t_*i*, δ_*i*) the observed time and censoring indicator. For DS1 (leukemia), the single covariate is the initial white blood cell count (p = 1). For DS3 (electronic components), no subject-level covariates are recorded in the original Aarset (1987) dataset; the PISN is therefore run in an intercept-only configuration (p = 0, bias term only) analogous to DS2. All three datasets are thus evaluated in the same framework, with the PISN reducing to a parametric GLD fit when covariates are absent.

The PISN maps x_*i* → (a_*i*, b_*i*, c_*i*) through

hℓ=ReLU(BN(Wℓhℓ-1+bℓ)), ℓ=1,2,(20)

(ai,bi,ci)T=Softplus(W3h2+b3),(21)

where BN(·) is batch normalization and Softplus enforces positivity of the GLD parameters. The individual survival prediction is Ŝ(t|xi)=exp(-ai/bi(ebit-1))(1+cit)-1. The architecture is illustrated in Figure 10.

The architecture uses two hidden layers of 64 units each (approximately 4,200 parameters for p = 1), which constitutes an overparameterized regime relative to DS1 (n = 33) and DS3 (n = 50). To mitigate overfitting, batch normalization is applied after every linear layer, an L₂ weight-decay penalty of 10⁻⁴ is included in the Adam optimizer, and all reported metrics are averages from 5-fold cross-validation. We acknowledge that a smaller architecture or stronger regularization may be preferable for very small samples (n < 50), and encourage practitioners to verify via cross-validated model selection. The results on these small datasets should be interpreted with appropriate caution: the observed improvements in C-index over parametric baselines (Table 6) are consistent but modest, and the confidence intervals are wide. The primary strength of the PISN on small datasets is the interpretability of the fitted GLD parameters remain meaningful, and the hazard regime diagnosis is directly readable from âb^≷ĉ2 rather than a claim of superior predictive performance in all small-sample settings.

Table 6. Survival prediction performance (5-fold cross-validation mean ± SD) of the PISN vs. benchmark survival models on DS1 (leukemia, n = 33), DS2 (COVID-19, n = 50), and DS3 [electronic component failures, n = 50; (Aarset, 1987)].

^†DS2 and DS3 contain no subject-level covariates; all models are run in intercept-only (marginal) configuration for these datasets. Bold values indicate the best-performing model for each metric (highest C-index, lowest IBS) within each dataset.

6.1.1 Computational cost

Training the PISN on DS1 (n = 33, p = 1) requires fewer than 60 s on a single CPU core (Python 3.10, PyTorch 2.1, Intel Core i7-1185G7). For DS3 (n = 50), the training time is under 90 s. Five-fold cross-validation adds a factor of five to these figures. The dominant cost is the Brent quantile solver called at each training step to verify the physics penalty; this can be parallelized across subjects in the future implementations. For comparison, GLD-MLE fitting via L-BFGS-B completes in under 2 s on the same hardware. The PISN is therefore computationally accessible for small-to-moderate clinical and engineering datasets on standard hardware. A larger benchmark dataset with subject-level covariates remains a priority for future work.

6.2 Physics-informed loss function

The total loss combines Equations 22–24. The total training loss is L=LNLL+γ1Lshape+γ2Lcal, where the censored negative log-likelihood is

LNLL=-1n∑i=1n[δilnĥ(ti|xi)+lnŜ(ti|xi)].(22)

The hazard shape penalty (from Theorem 0.0.4) is

Lshape=1n∑i=1nmax(0,σκ(âib^i-ĉi2)·h′(ti;âi,b^i,ĉi))2,(23)

where σ_κ(u) = tanh(κu) with κ = 20 is a smooth, everywhere-differentiable approximation to sign(u). This ensures that Lshape and its gradient with respect to (âi,b^i,ĉi) are well-defined at every point in the parameter space, consistent with the differentiable-penalty claim of Section 6. For κ = 20, σ_κ(u) deviates from sign(u) by less than 10⁻⁸ for |u|>1. The penalty penalizes h′(ti) whenever its sign contradicts the predicted IFR or bathtub regime.

The population calibration loss is

Lcal=|T|-1∑t∈T(ŜKM(t)-Ŝ¯(t))2.(24)

The PISN is implemented in PyTorch (Paszke et al., 2019), optimized with Adam (Kingma and Ba, 2015). Hyperparameter details are provided in Section 6.3.

6.3 Hyperparameter selection and ablation

The penalty weights γ₁ and γ₂ are selected by 5-fold cross-validation on the training split over the grid γ₁, γ₂∈{0.01, 0.1, 1.0}, choosing the combination that minimizes the mean validation concordance index (C-index). Two hidden layers with 64 units each were sufficient for all datasets considered here.

To assess the contribution of each penalty term, Table 7 reports an ablation study comparing four PISN configurations on DS1 and DS3, the two datasets with the most distinct hazards regimes. Removing the shape penalty (γ₁ = 0) leads to a consistent reduction in C-index and an increase in IBS, indicating that the hazard-shape constraint from Theorem 0.0.4 contributes meaningfully to predictive performance. The full PISN with both penalties active yields consistently improved performance across all configurations examined.

Table 7. Ablation study: C-index (higher is better) and IBS (lower is better) for four PISN configurations on DS1 (leukemia) and DS3 (electronic components).

Bold values indicate the best-performing model for each metric (highest C-index, lowest IBS) within each dataset.

6.4 PISN benchmark results

Table 6 compares the PISN against standard parametric and deep learning survival baselines on all three datasets: leukemia, COVID-19, and electronic component failures. All metrics are reported as mean ± standard deviation across 5 folds. The hazard regime for each dataset is identified by the fitted condition âb^≷ĉ2 from Theorem 0.0.4. Since DS2 records only marginal survival times (no subject-level covariates), the PISN, and all covariate-driven baselines are run in their covariate-free (intercept-only) configuration for DS2; the GLD-MLE entry therefore coincides with the parametric fit reported in Section 7 and serves as the natural reference.

The PISN yields consistently higher C-index and lower IBS values across all three datasets. The physics penalty correctly diagnoses the increasing hazard (IFR) regime on DS1, a near-boundary decreasing hazard (DFR) regime on DS2, consistent with the early excess mortality pattern in acute COVID-19 hospitalization (Ahamed et al., 2023), and the bathtub (UBT) regime on DS3, all consistent with Theorem 0.0.4 and the information-criteria analysis of Section 7. The RSF baseline performs comparably to the deep learning models, confirming the competitiveness of non-parametric approaches on small tabular datasets.

7 Applications to real data

7.1 Datasets and goodness-of-fit criteria

Three datasets from recently published studies are analyzed. All information criteria and the Kolmogorov-Smirnov statistic are reported:

AIC=-2ℓ^+2k, AICC=AIC+2k(k+1)n-k-1,

BIC=-2ℓ^+kln n, D=supx|Fn(x)-F^(x)|,(25)

All four criteria (Equation 25) are reported in Tables 1–3. where k is the number of estimated parameters.

Dataset 1 (DS1): Survival times (weeks) of n = 33 patients with acute myelogenous leukemia from Feigl and Zelen (1965), recently used as a benchmark by (Eghwerido et al. (2021) and (Ihtisham et al. 2019). Descriptive statistics: mean = 40.88 wk, SD = 45.99 wk, skewness = 1.69, kurtosis = 4.87. The covariate used in the PISN analysis is the log-transformed initial white blood cell count recorded in the original study.

Dataset 2 (DS2): Survival days of n = 50 hospital COVID-19 patients from the supplementary data of (Ahamed et al. 2023). A publicly mirrored version of the same data is available at https://www.kaggle.com/datasets/meirnizri/covid19-dataset. All analyses use the data as reported in (Ahamed et al. 2023). Descriptive statistics: mean = 42.92 days, SD = 55.80 days, skewness = 1.67, kurtosis = 4.59. No subject-level covariates are available; the PISN is run in intercept-only configuration.

Dataset 3 (DS3): Failure times (hours) of n = 50 electronic components from the accelerated life test reported by Aarset (1987). This is the canonical benchmark for bathtub-shaped hazard analysis; the bathtub empirical hazard rate was established in the original paper and confirmed in many subsequent reliability studies. No subject-level covariates are recorded; the PISN is run in intercept-only configuration. Descriptive statistics: mean = 45.69 h, SD = 32.51 h, skewness = −0.14, kurtosis = 1.88.

The competing distributions are listed in Table 8.

Table 8. Competing distributions used for goodness-of-fit comparison. All parameter estimates are obtained by MLE.

7.2 Diagnostic plots

Figure 11 shows the fitted GLD density overlaid on histograms for each dataset. Figure 12 presents Q-Q plots in which points close to the 45° line confirm adequate fit. Figure 13 shows observed vs. expected frequency bar charts, providing a discrete goodness-of-fit visual for each dataset. All three datasets show satisfactory agreement between the observed data and the fitted GLD.

8 Conclusion

The Gompertz-Lomax Distribution establishes a new three-parameter lifetime model whose additive composite hazard h(x) = ae^*bx*+c/(1+cx) provides, to the best of the authors' knowledge, the first analytically tractable combination of a Gompertz IFR component with a Lomax DFR component yielding a fully closed-form survival function. The hazard shape theorem establishes the single condition ab≷c² that separates IFR from bathtub behavior, and this same condition serves as the differentiable physics constraint in the PISN loss function, implemented via a smooth tanh approximation to ensure end-to-end gradient propagation. The empirical results confirm that the GLD achieves among the lowest K-S statistics and the most favorable information criteria on all three benchmark datasets, with a decisive margin over all competing models on the Aarset electronic component dataset, the only benchmark whose known bathtub hazard structure the GLD can analytically represent. On the leukemia and COVID-19 datasets, the GLD leads on every information criterion, confirming that the composite hazard structure captures patterns missed by simpler two-parameter models.

Despite these contributions, several limitations should be acknowledged. First, while S(x), f(x), and h(x) are fully closed-form, several downstream quantities, notably the quantile function, moments, mean residual life, and stress-strength reliability require numerical evaluation. This is explicitly noted in Section 3 and does not diminish the distributional tractability, but should be considered in time-critical applications. Second, the three-parameter structure increases estimation sensitivity in small samples (n < 30), as reflected by the higher MSE values at n = 25 in Table 4; practitioners working with very small datasets may prefer to fix one parameter or use Bayesian estimation with informative priors. Third, the PISN architecture assumes that the GLD parametric form is appropriate for the underlying survival mechanism; misspecification of this assumption will propagate into biased parameter estimates, and a goodness-of-fit test for the GLD family should be conducted before deployment. Fourth, the PISN evaluation datasets (DS1: n = 33; DS2, DS3: n = 50) are small. The observed performance improvements over parametric baselines (Table 6) are consistent but modest; generalizability to larger datasets with richer covariates structures cannot be established from the present evaluation alone, and replication on larger cohorts is strongly encouraged. Fifth, the stress-strength reliability integral does not reduce to a simple closed form when the two populations have different shape parameters b₁≠b₂, limiting the direct analytical treatability of that setting.

Future directions include a zero-inflated GLD for excess-zero failure count data, a bivariate GLD via the Marshall-Olkin copula for competing-risk settings, Bayesian estimation using Hamiltonian Monte Carlo for small clinical samples, and time-varying covariate extensions of the PISN through recurrent or attention-based encoders. Evaluation of larger survival benchmarks such as SEER or MIMIC remain a priority.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Qun Dai, Nanjing University of Aeronautics and Astronautics, China

Reviewed by: Abdulrahman Alshammari, Jouf University, Saudi Arabia

Akbar Khan, Government Post Graduate College Kohat, Pakistan