A novel Weighted Exponential-Lindley distribution with Bayesian estimation under progressive Type-II censoring and a Bayesian-Calibrated Variational Autoencoder for survival prediction

1 Introduction

Precise probabilistic modeling of lifetime data occupies a central position in reliability engineering, clinical survival analysis, and risk management (Lawless, 2011). Traditional one-parameter lifetime distributions such as the exponential and Lindley distributions are valued for their simplicity, but are structurally restricted in the range of dispersion and tail behavior they can represent, since a single parameter must simultaneously govern both location and shape: fixing θ fixes the mean, variance, coefficient of variation, and hazard level all at once, leaving no free parameter with which to adjust the distribution's shape independently of its scale, once a dataset's central tendency has pinned down θ.

1.1 The critical role of weighted distributions

Weighted distributions form a foundational and still-active branch of distribution theory that directly motivates the present work. Rao's original formulation (Rao, 1965) showed that many observed sampling schemes do not draw observations with uniform probability from an underlying population, but rather record each unit with probability proportional to a weight function w(x), so that the observed density f^*(x) ∝ w(x)f(x) can differ systematically from the population density f(x); length-biased sampling, w(x) = x, is the special case most relevant to lifetime data, since units with longer lifetimes are mechanically more likely to be sampled under many observational schemes. Patil and Rao (1978) formalized this framework and demonstrated its relevance across wildlife population and human family-size studies, and it has since become one of the standard mechanisms for constructing new, more flexible lifetime distributions from simple parametric bases: applying a weight function to an existing density adds shape flexibility while preserving a direct probabilistic interpretation of the added parameter, in contrast to purely mathematical generator families that add parameters without such an interpretation. The weighted exponential distribution of Gupta and Kundu (2009), obtained by applying Azzalini's skewing construction to the exponential density is a widely cited example of this mechanism in the reliability literature, and the approach remains active and productive today, with recent weighted extensions proposed across a wide range of base distributions, including the weighted exponential-Gompertz distribution (Abd El-Bar and Ragab, 2019), the weighted log-Kumaraswamy distribution (Ragab et al., 2026), and weighted extensions of the Ailamujia distribution for radiation data (Alballa et al., 2026). A closely related route to added shape flexibility, the power transformation Y = T^1/δ applied to a base variable T, has likewise been used recently to construct flexible two-parameter lifetime models such as the power Obulezi distribution (Zhao et al., 2026); a further route combines this power transformation with domain truncation, as in the truncated power Chris–Jerry distribution (Hassan et al., 2025), which restricts a power-transformed base model to a bounded interval to better accommodate proportion- or rate-type data. While mechanistically distinct from weighting, both routes pursue the same broad goal of extending a simple base distribution's shape and hazard repertoire with one or two additional parameters, and we return to this family, alongside the weighted comparators above, when motivating the WEL construction below. The WEL distribution developed in this article belongs squarely to this tradition: it is constructed (Section 2.1) by applying a length-biased weight w(x) = x^*c* to an exponential-Lindley mixture base, and Section 7 compares it not only against unweighted competitors, but against multiple other members of the weighted-distribution class, including the weighted Lindley distribution of Ghitany et al. (2011) and the weighted exponential distribution of Gupta and Kundu (2009).

Motivation The specific motivation for this article follows directly from the discussion above: existing weighted lifetime distributions, including the weighted Lindley distribution, apply their weight function to a single base component, so the weight parameter and the scale parameter necessarily interact within that one component and cannot be tuned to independent statistical targets (made precise below). We construct the WEL distribution to break this coupling by first forming a two-component exponential-Lindley mixture that shares a single scale parameter, and only then applying a length-biased weight across the mixture as a whole, so that the weight parameter's effect on dispersion and tail behavior is structurally separated from the scale parameter's role as the asymptotic hazard level. This construction, its consequences for hazard-rate shape, and its practical value relative to the weighted Lindley and other weighted distributions on real data are the subject of the remainder of this article.

The Lindley distribution, introduced by Lindley (1958) as a counterexample in fiducial statistics has generated an extensive family of one- and two-parameter extensions due to its interpretable mixture representation and closed-form survival function (Ghitany et al., 2008; Shanker, 2015, 2016, 2017; Shukla, 2018), and this generalization effort remains active: recent extensions include Marshall-Olkin and generalized-Lindley constructions (Algarni, 2021; Abouammoh et al., 2015), mixture-based extensions combining Lindley and gamma or exponential components in the same spirit as the WEL construction below (Elgarhy et al., 2024), heavy-tailed power-transformed extensions aimed at actuarial applications (Elgarhy et al., 2025), and broader generalized-G families built by compounding generator families (Gabanakgosi and Otlaadisa, 2024). Among the Lindley extensions, the two-parameter weighted Lindley distribution proposed by Ghitany et al. (2011) introduced a length-biased weighting mechanism into this family, demonstrating that the resulting model fits several benchmark survival datasets better than the unweighted Lindley and exponential distributions. This is the specific gap the present article targets: in the weighted Lindley density, the weight parameter α and the scale parameter θ enter the same single Lindley mixture component multiplicatively, so that adjusting α to fit a dataset's dispersion necessarily perturbs the implied hazard level governed by θ, and vice versa; the two parameters cannot be tuned to independently target dispersion/tail weight and hazard level. The exponential-Lindley base mixture we construct in Section 2.1, by contrast, retains a single shared scale θ across both mixture components by design, so that the length-biased weight c we subsequently apply (Section 2.1) modifies dispersion and tail behavior while θ remains cleanly interpretable as the asymptotic hazard level (Theorem 2.2); this structural separation, not merely the addition of a third parameter, is the specific limitation of the exponential-Lindley family that motivates the WEL construction.

A complementary limitation characterizes contemporary deep learning survival models, which have grown rapidly enough that several recent surveys now track the field systematically (Wiegrebe et al., 2024). DeepSurv (Katzman et al., 2018) extends the Cox proportional hazards model with a neural network but retains the proportional-hazards assumption and provides no posterior uncertainty over the individual predictions. DeepHit (Lee et al., 2018) discretises the time axis and produces point estimates of the survival function without an associated probability distribution over those estimates. Deep Survival Machines (Nagpal et al., 2021) is closely related in spirit to the present work: it also parametrises individual survival distributions (mixtures of Weibull or Log-Normal components) via a neural network, but its component distributions are fitted by direct maximum likelihood rather than full variational inference with an explicit latent space and posterior predictive uncertainty. The Survival Variational Autoencoder of Apellániz et al. (2024) (SAVAE) and the more recent DySurv model of Mesinovic et al. (2024), which uses a conditional variational autoencoder for dynamic survival prediction, are the closest prior work combining variational inference with parametric or semi-parametric survival modeling; SAVAE's authors employ an exponential or Weibull parametric family as the decoder output and explicitly identify the construction of calibrated uncertainty intervals for the individual predictions as an open direction for future work. More broadly, existing parametric survival models and deep learning survival models occupy largely separate methodological territories: the former offer closed-form interpretability but limited flexibility, while the latter offer empirical flexibility but typically sacrifice both interpretability and principled uncertainty quantification.

The present article makes two contributions, and we state their novelty explicitly. The first is the Weighted Exponential-Lindley (WEL) distribution, a three-parameter lifetime model constructed by applying a length-biased weight function w(x) = x^*c* to a mixture of an exponential component and a Lindley component sharing a common scale parameter θ—a shared-scale construction that, to our knowledge, has not been used to decouple the hazard level from dispersion control within the exponential-Lindley family, specifically, as distinct from the many existing generalized-Lindley constructions cited above, which adds shape flexibility without this particular separation property. We show that this construction yields a two-component gamma mixture representation, from which we prove, using Glaser's failure-rate criterion applied directly to the WEL density (Section 2.5) rather than by appeal to an external mixture theorem, that the WEL hazard rate is strictly increasing on its entire parameter space, converging to θ as x → ∞. The weight parameter c does not alter this qualitative IFR shape, but we show that it provides a continuous tunable, independent control over the coefficient of variation, skewness and kurtosis of the distribution, a property we verify both analytically and through real-data model fits in which the additional parameter measurably improves the maximized log-likelihood relative to the unweighted Lindley distribution, and, on some but not all real datasets examined, relative to the weighted Lindley distribution as well (Section 7).

The second contribution is the Bayesian-Calibrated Variational Autoencoder Survival Model (BC-VAESM), a deep generative architecture whose decoder output layer produces individual, covariate-dependent WEL parameters rather than a single fixed population-level parameter vector combining, within a single architecture, the full variational-inference latent-space structure of SAVAE/DySurv with a flexible IFR parametric family in place of the exponential or Weibull decoders used in that prior work. We train the BC-VAESM end-to-end via a censored negative log-likelihood combined with a Kullback-Leibler regularization term, implement the WEL log-density and log-survival functions in a numerically stable, fully differentiable form suitable for back-propagation, and evaluate the resulting model through five-fold cross-validation on real clinical data. To isolate the genuine contribution of the WEL parametric family within the deep learning architecture, we additionally train an architecturally identical variational autoencoder whose decoder instead outputs a two-parameter Weibull distribution, holding the encoder, latent dimension, training procedure, and all hyperparameters fixed.

In addition, we develop a Bayesian estimation of the WEL distribution under progressive Type-II censoring, a censoring scheme of established relevance in clinical trials and accelerated life testing, where units are removed from the study at intermediate failure times (Balakrishnan and Aggarwala, 2000). Bayesian estimators under squared error, LINEX, and general entropy loss functions are derived and compared against maximum likelihood and method-of-moments estimators through a Monte Carlo simulation study spanning three parameters configurations and sample sizes from n = 25 to n = 300 (Section 5), with explicit reporting of the estimator convergence rates, confidence/credible interval coverage, and per-replication computation time.

The remainder of this article is organized as follows. Section 2 defines the WEL distribution, establishes its mixture representation, and proves the hazard rate theorem, including an explicit statement of the regularity conditions and an identifiability result (Section 4.1.1). Section 3 derives the complete set of distributional properties. Section 4 develops maximum likelihood, method-of-moments, and Bayesian estimation under progressive Type-II censoring, with the score functions and observed information matrix given in Appendix 1. Section 5 reports the Monte Carlo simulation study. Section 6 introduces the BC-VAESM architecture and training procedure. Section 7 applies both the WEL distribution and the BC-VAESM to four real datasets, with full transparency regarding cases in which the WEL distribution does not outperform simpler competitors. Section 9 concludes with a discussion of limitations and future directions.

2 The Weighted Exponential-Lindley distribution

2.1 Motivation and construction

Let f₀(x; θ, α) denote the density of the Exponential-Lindley mixture:

f0(x;θ,α)=αθe-θx+(1-α)θ2θ+1(1+x)e-θx,x>0,θ>0,α∈(0,1).(1)

The first component is the exponential density with scale θ; the second is the Lindley density (Lindley, 1958) with parameter θ. The mixing weight α controls the relative contribution of the two components.

Following the length-biased weighting framework of Patil and Ord (1976), we apply the weight function w(x) = x^*c*, c > 0, to obtain the weighted density

f(x;θ,α,c)=xcf0(x;θ,α)μc(θ,α,c), x>0,(2)

where μc=𝔼f0[Xc] is the normalizing constant. Writing the numerator of Equation 2 as x^*c*e^−θ*x*A + x^{*c* + 1}e^−θ*x*B with

A=αθ+(1-α)θ2θ+1, B=(1-α)θ2θ+1,(3)

and integrating term by term using ∫0∞xke-θxdx=Γ(k+1)/θk+1 gives the normalizing constant in closed form:

μc(θ,α,c)=AΓ(c+1)θc+1+BΓ(c+2)θc+2.(4)

Definition 2.1 (WEL Distribution). *A continuous random variable X follows the WEL distribution with parameter vector *ψ* = (θ, α, c), written X ~ WEL(θ, α, c), if its probability density function is*

f(x;θ,α,c)=xce-θxμc[αθ+(1-α)θ2(1+x)θ+1],x>0,θ>0,α∈(0,1),c>0,(5)

where μ_c is given by Equation 4.

2.2 Mixture representation

Proposition 2.1 (Mixture Form). The WEL density Equation 5 can be written as the two-component gamma mixture

f(x;θ,α,c)=p1g(x;n1,θ)+p2g(x;n2,θ),(6)

*where g(x; n, θ) denotes the Gamma(n, θ) density with shape n and rate θ, the shape parameters are n₁ = c + 1 and n₂ = c + 2, and the mixing weights are*

p1=AΓ(n1)/θn1μc, p2=BΓ(n2)/θn2μc,(7)

with A, B as in Equation 3, satisfying p₁ + p₂ = 1, p₁, p₂ > 0.

Proof. The numerator of Equation 5 is, by direct expansion of the bracket, x^*c*e^−θ*x*A + x^*c*+1e^−θ*x*B. Writing each term as a constant multiple of the corresponding Gamma(n_*j*, θ) density, g(x;nj,θ)=θnjxnj-1e-θx/Γ(nj), and normalizing by μ_*c*, which by construction equals the sum of the two unnormalised integrals (Equation 4), yields Equations 6–7 with p₁ + p₂ = 1 following immediately since μ_*c* is exactly the sum of the two numerators evaluated at their respective normalizing constants. □

We verified Equations 6–7 numerically: for every parameter combination examined, p₁ + p₂ equals 1 to machine precision, and the implied cumulative distribution function agrees with direct numerical integration of Equation 5 to at least five decimal places across the support.

2.3 Special cases

• As c → 0⁺: μ_*c* → 1 and f(x; θ, α, 0) = f₀(x; θ, α), the base Exponential-Lindley mixture Equation 1.
• As α → 1, c → 0⁺: f → θe^−θ*x*, the Exponential distribution.
• As α → 0, c → 0⁺: f → [θ²/(θ + 1)](1 + x)e^−θ*x*, the Lindley distribution (Lindley, 1958).

2.4 Survival and hazard functions

Using Proposition 2.1, the survival function is

F̄(x;ψ)=p1Q(n1,θx)+p2Q(n2,θx),(8)

where Q(n, z) = Γ(n, z)/Γ(n) is the regularized upper incomplete gamma function, and the cumulative distribution function is F(x;ψ)=1-F̄(x;ψ). The hazard rate function is

h(x;ψ)=f(x;ψ)F̄(x;ψ)=xce-θx[A+Bx]/μcp1Q(n1,θx)+p2Q(n2,θx).(9)

2.5 Hazard shape theorem

Theorem 2.2 (Hazard Shape Characterization). Let X ~ *WEL(θ, α, c) with θ* > 0, α ∈ (0, 1), c > 0. Then the hazard rate function h(x; *ψ*) is strictly increasing on (0, ∞), with

h(0+;ψ)=0 and limx→∞h(x;ψ)=θ.(10)

That is, the WEL distribution belongs to the increasing-failure-rate (IFR) class for every admissible parameter value, with the scale parameter θ acting as the asymptotic (long-run) hazard level.

Proof. We use Glaser's failure-rate criterion (Glaser, 1980), which applies directly to the WEL density and does not require invoking an external mixture-monotonicity result. Define η(x) = −f′(x; *ψ)/f(x*; *ψ). Glaser's theorem states that if η′(x*) > 0 for all x > 0 and f(0⁺) is finite (or zero), then the hazard rate h(x; *ψ) is strictly increasing on (0, ∞). From Equation 5, dropping the normalizing constant μ_c* (which does not depend on x and hence does not affect η),

ln f(x;ψ)=const+cln x-θx+ln(A+Bx),(11)

with A, B as in Equation 3. Differentiating,

f′(x;ψ)f(x;ψ)=cx-θ+BA+Bx, so η(x)=θ-cx-BA+Bx.(12)

Differentiating again,

η′(x)=cx2+B2(A+Bx)2.(13)

Since c > 0 by assumption, the first term of Equation 13 is strictly positive for every x > 0, and the second term is non-negative (strictly positive whenever α < 1, i.e. whenever the Lindley component has positive weight, since then B > 0); hence η′(x) > 0 for all x > 0 regardless of the value of α ∈ (0, 1). By Proposition 2.1, f(x; *ψ*) = p₁g(x; n₁, θ) + p₂g(x; n₂, θ) with n₁ = c + 1 > 1 and n₂ = c + 2 > 1 whenever c > 0, so g(0;n_*j*, θ) = 0 for j = 1, 2 and therefore f(0⁺; *ψ*) = 0, satisfying Glaser's boundary condition. By Glaser's theorem, h(x; *ψ*) is therefore strictly increasing on (0, ∞) for every θ > 0, α ∈ (0, 1), c > 0—i.e. for every point in the admissible parameter space, with no additional restriction on α or θ beyond their defining ranges. Since f(0⁺; *ψ*) = 0, h(0+;ψ)=f(0+;ψ)/F̄(0;ψ)=0. For the limit as x → ∞, from Equation 9, h(x;ψ)=xce-θx(A+Bx)/[μcF̄(x;ψ)]; applying L'Hôpital's rule to the ratio f(x)/F̄(x) as x → ∞ (both → 0), or equivalently, invoking the standard asymptotic hazard result for gamma-type tails limx→∞g(x;n,θ)/Ḡ(x;n,θ)=θ (which holds for both mixture components since they share the rate θ), together with the fact that a mixture of two hazards sharing the same limit θ itself converges to that common limit, gives limx→∞h(x;ψ)=θ. □

Remark 2.1 (Relation to the Block–Savits IFRA closure result). *An earlier version of this proof invoked the mixture-monotonicity result of Block and Savits (1976) directly; that result concerns closure of the weaker IFRA (increasing failure rate average) class under mixtures and does not by itself establish the strict pointwise monotonicity of the hazard rate (as opposed to its average) claimed in Theorem 2.2. We have therefore replaced that argument with the self-contained Glaser-criterion proof above, which establishes the stronger IFR property directly from the WEL density and makes explicit exactly which assumptions* (c > 0, θ > 0, α ∈ (0, 1)) are used at each step.

We verified Theorem 2.2 numerically by evaluating h(x; *ψ*) over a fine grid in x for dozens of parameter combinations spanning four orders of magnitude in θ, α, and c: in every case examined, h(x) was strictly increasing and converged numerically to θ. Figure 1 illustrates this behavior for two representative slices of the parameter space, one varying θ and one varying c; every curve rises monotonically toward its dashed asymptote, exactly as Theorem 2.2 predicts, regardless of which parameter is varied.

Remark 2.2. *Theorem 2.2 shows that the WEL distribution belongs to the IFR class throughout its parameter space, with the scale parameter θ interpretable as the limiting hazard level a property not shared by the one-parameter Lindley distribution, whose hazard rate is instead bounded above by θ²/(θ + 1) < θ. While the WEL distribution therefore does not admit bathtub or upside-down-bathtub hazard shapes, its principal contribution lies in the weight parameter c, which, as demonstrated in Section 3.2, provides an independent and continuously tunable control over the coefficient of variation, skewness, and kurtosis at fixed θ, substantially extending the range of dispersion and shape behavior attainable within the Exponential-Lindley mixture family. This dispersion-control role of c, combined with the closed-form normalizing constant Equation 4 that renders exact likelihood-based and Bayesian inference tractable even under censoring (Section 4), constitutes the principal theoretical contribution of the WEL family. We return to this point in Section 7, where the weight parameter is shown to yield a measurable improvement in log-likelihood over the unweighted Lindley and weighted Lindley models on real data, even though it does not change the qualitative hazard shape.*

2.6 Probability density function: graphical illustration

Figure 2 displays the WEL density for a range of θ and c values, illustrating the right-skewed, unimodal shape typical of the family and the manner in which c governs both the location of the mode and the heaviness of the right tail.

3 Statistical properties

3.1 Moments and generating functions

Using the mixture representation (Equation 6) and the gamma raw-moment formula 𝔼[Gamma(n, θ)^*r*] = Γ(n + r)/[θ^*r*Γ(n)], the r-th raw moment of the WEL distribution is

μr′=𝔼[Xr]=1μc[AΓ(c+r+1)θc+r+1+BΓ(c+r+2)θc+r+2],r=1,2,3,…(14)

with A, B as in Equation 3 and μ_*c* as in Equation 4. The variance is σ2=μ2′-(μ1′)2. We verified Equation 14 numerically against direct numerical integration of x^*r*f(x) for r = 1, 2, 3, 4 across the same parameter grid used to verify the mixture representation, finding agreement to at least five significant figures throughout.

The moment generating function follows from Equation 6 and the gamma MGF 𝔼[e^*tX*] = (θ/(θ − t))^*n* for Gamma(n, θ):

MX(t)=p1(θθ-t)n1+p2(θθ-t)n2, t<θ,(15)

and the characteristic function is φ_*X*(t) = M_*X*(it).

3.2 Skewness, kurtosis, coefficient of variation, and dispersion index

The standard shape measures follow from the central moments μr=𝔼[(X-μ1′)r], obtained from the raw moments (Equation 14) via the binomial expansion μr=∑k=0r(rk)(-μ1′)r-kμk′:

β1=μ3μ23/2, β2=μ4μ22, CV=μ2μ1′, γ=μ2μ1′.(16)

Figure 3 displays these measures as functions of c at fixed θ = 1.0, α = 0.5, computed directly from Equation 14. All four measures decrease monotonically and smoothly as c increases, confirming that c acts as an independent dispersion-control parameter: larger c produces a distribution that is less variable, less right-skewed, and less heavy-tailed relative to its mean, approaching Gaussian-like shape measures as c grows, while θ separately controls scale.

3.3 Mean residual life

The mean residual life (MRL) function, which measures the expected remaining lifetime given survival to age x₀, is

m(x0)=1F̄(x0)∫x0∞F̄(t)dt,(17)

with m(0)=μ1′ by definition. Section 3.4 below derives the general s-th conditional moment of the WEL distribution and shows that m(x₀) admits the closed form m(x0)=φ1(x0)/F̄(x0;ψ)-x0 (Equations 19–21), so the numerical integration of Equation 17 used to produce Figure 4 below serves only as an independent numerical check on that closed-form expression, not as the primary means of evaluating it. Consistent with the IFR property established in Theorem 2.2, m(x₀) is decreasing in x₀ in every configuration examined, a well-known consequence of the IFR property (Shaked and Shanthikumar, 2007). Figure 4 confirms this monotone decrease across the parameter configurations shown.

3.4 Conditional moments

Conditional moments generalize the mean residual life of Section 3.3 to arbitrary order and to both tails of the distribution, and they underlie several of the reliability and inequality measures used elsewhere in this article (the mean residual life itself, the mean inactivity time introduced below, and the Bonferroni and Lorenz curves of Section 3.10). We therefore derive them explicitly here, rather than leaving their construction implicit.

Because the mixture representation of Proposition 2.1 writes f(x; *ψ*) as a two-component gamma density, any truncated integral of x^*s*f(x; *ψ) reduces to a weighted sum of truncated gamma moments, each available in closed form through the (unregularised) incomplete gamma function. Specifically, for a single Gamma(n*, θ) component,

∫t∞xsg(x;n,θ)dx=θ-sΓ(n+s,θt)Γ(n),∫0txsg(x;n,θ)dx=θ-sγ(n+s,θt)Γ(n),(18)

where Γ(·, ·) and γ(·, ·) denote the upper and lower incomplete gamma functions, respectively (so that Γ(n + s, θt) + γ(n + s, θt) = Γ(n + s) for every t ≥ 0). Applying Equation 18 to each of the two components in Equation 6, with n₁ = c + 1 and n₂ = c + 2 as before, gives the s-th upper and lower partial moments of X,

φs(t)=∫t∞xsf(x;ψ)dx =θ-s[p1Γ(n1+s,θt)Γ(n1)+p2Γ(n2+s,θt)Γ(n2)],(19)

χs(t)=∫0txsf(x;ψ)dx =θ-s[p1γ(n1+s,θt)Γ(n1)+p2γ(n2+s,θt)Γ(n2)],(20)

for s = 1, 2, 3, …, with φs(t)+χs(t)=μs′ (the raw moment of Equation 14) for every t ≥ 0, as required by construction, since the two truncated regions partition the entire support. The s-th conditional moments follow immediately by dividing by the corresponding tail probability:

𝔼[Xs∣X>t]=φs(t)F̄(t;ψ), 𝔼[Xs∣X≤t]=χs(t)F(t;ψ),(21)

with F̄(t;ψ) as in Equation 8. Two special cases recover quantities used elsewhere in the article. Setting s = 1 in the first expression of Equation 21 and subtracting t reproduces the mean residual life m(t)=φ1(t)/F̄(t;ψ)-t of Equation 17, now available in closed form rather than requiring numerical integration; and the mean inactivity time (the expected elapsed time since the preceding renewal, given that it has already occurred by time t),

ι(t)=t-𝔼[X∣X≤t]=t-χ1(t)F(t;ψ), t>0,(22)

follows directly by setting s = 1 in the second expression of Equation 21. We verified Equations 19–22 numerically against direct numerical integration of x^*s*f(x; *ψ) over [t*, ∞) and [0, t] for s = 1, 2 and a grid of t and parameter values, finding agreement to at least five significant figures throughout, consistent with the level of agreement already reported for the raw moments in Section 3.1.

3.5 Quantile function

The p-th quantile Q(p) satisfies F̄(Q(p))=1-p; there is no closed-form inverse, so Q(p) is computed numerically via root finding with the exponential quantile −log(1 − p)/θ as the starting value. The median is M = Q(0.5).

3.6 Order statistics

Let X₍₁₎ ≤ ⋯ ≤ X_(*n*) be the order statistics of an i.i.d. sample of size n from the WEL distribution. The density of the k-th order statistic is

fX(k)(x)=n!(k-1)!(n-k)![F(x)]k-1[F̄(x)]n-kf(x),(23)

which, upon substituting Equations 5 and 8 and expanding [F(x)]^*k*−1 via the binomial theorem, yields a finite series of WEL survival and density terms.

3.7 Stochastic ordering

Theorem 3.1 (Stochastic Ordering). Let X_*i* ~ WEL(θ_*i*, α, c), i = 1, 2, with θ₁ ≥ θ₂. Then X₁ ≤_lr X₂, implying X₁ ≤_hr X₂, X₁ ≤_mrl X₂, and X₁ ≤_st X₂.

Proof. The likelihood ratio is

f1(x)f2(x)=μc,2μc,1·αθ1+(1-α)θ12(1+x)/(θ1+1)αθ2+(1-α)θ22(1+x)/(θ2+1)·e-(θ1-θ2)x.

For θ₁ > θ₂, differentiating the logarithm of this ratio with respect to x yields a strictly negative expression (the derivative of the log of the rational factor is bounded while the linear term −(θ₁ − θ₂)x dominates), establishing X₁ ≤_lr X₂. The implication chain ≤_lr⇒ ≤_hr⇒ ≤_mrl⇒ ≤_st follows from standard stochastic order theory (Shaked and Shanthikumar, 2007). □

3.8 Stress-strength reliability

Let X ~ WEL(θ₁, α₁, c) (strength) and Y ~ WEL(θ₂, α₂, c) (stress) be independent. The stress-strength reliability is

R=P(Y<X)=∫0∞FY(x)fX(x)dx,(24)

which, substituting the mixture representations of F_*Y* and f_*X* from Proposition 2.1, reduces to a finite sum of terms of the form ∫0∞[1-Q(n2,j,θ2x)]g(x;n1,i,θ1)dx, each evaluable via the regularized incomplete beta function after the substitution u = θ₁x/(θ₁ + θ₂).

3.9 Entropy measures

Entropy quantifies the uncertainty, or average information content, carried by a random variable, and is a standard descriptive tool in the lifetime-distribution literature alongside the moment-based shape measures of Section 3.2. We derive the Rényi entropy of the WEL distribution in full below, together with its Shannon and Tsallis counterparts; the latter two follow almost immediately from the same intermediate integral and are therefore obtained at essentially no extra analytical cost once the Rényi derivation is in hand.

Rényi entropy. For δ > 0, δ ≠ 1, the Rényi entropy (Rényi, 1961) is defined as

TR(δ)=11-δln[∫0∞f(x;ψ)δdx].(25)

To evaluate the integral in Equation 25, raise the WEL density (Equation 5) to the power δ,

f(x;ψ)δ=μc-δxδce-δθx(A+Bx)δ,(26)

and expand the polynomial factor (A + Bx)^δ using the generalized binomial series (A+Bx)δ=∑k=0∞(δk)Aδ-kBkxk. Substituting this expansion into Equation 26 and integrating term by term against the gamma-type kernel, using ∫0∞xδc+ke-δθxdx=Γ(δc+k+1)/(δθ)δc+k+1, gives

∫0∞f(x;ψ)δdx=1μcδ∑k=0∞(δk)Aδ-kBkΓ(δc+k+1)(δθ)δc+k+1,(27)

which, substituted back into Equation 25, yields the closed-form series representation

TR(δ)=11-δln[1μcδ∑k=0∞(δk)Aδ-kBkΓ(δc+k+1)(δθ)δc+k+1].(28)

The series in Equation 28 converges rapidly because successive binomial coefficients decay once the summation index k exceeds δ; truncating at 40 terms (used throughout this article) reproduces the value obtained by direct numerical integration of Equation 25 to at least four significant figures over the parameter grid used elsewhere in the article.

Shannon entropy The Shannon entropy, H(X) = −𝔼[ln f(X; *ψ*)], arises as the limiting case limδ→1TR(δ) of the Rényi family (Rényi, 1961), and can equivalently be written directly from Equation 5 as

H(X)=lnμc-c𝔼[ln X]+θμ1′-𝔼[ln(A+BX)],(29)

where μ1′=𝔼[X] is available in closed form from Equation 14, while 𝔼[ln X] and 𝔼[ln(A + BX)] are evaluated by one-dimensional numerical integration against Equation 5, since neither admits an elementary closed form for general A, B, c. A larger value of H(X) indicates a flatter, less predictable density, and a smaller value indicates a more concentrated one, so H(X) provides a single-number summary that corroborates the dispersion trends already visible in the coefficient of variation and kurtosis of Section 3.2.

Tsallis entropy The Tsallis entropy (Tsallis, 1988), a non-extensive generalization of Shannon entropy, is widely reported alongside Rényi's measure in the weighted-distribution literature, is defined for δ > 0, δ ≠ 1, as

ST(δ)=1δ-1[1-∫0∞f(x;ψ)δdx],(30)

and follows immediately from the integral Equation 27 already derived for the Rényi entropy, without further computation. Since S_*T*(δ) and T_*R*(δ) are both monotonic transformations of the same underlying integral (Equation 27), they rank parameter configurations by dispersion identically, differing only in the functional form used to rescale that integral into an entropy value.

3.10 Bonferroni and Lorenz curves

The Lorenz and Bonferroni curves, which characterize the lifetime inequality, are

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

where ϕ1(t0)=∫0t0xf(x)dx is the first incomplete moment, computable from the mixture representation using the lower incomplete gamma function. The Gini coefficient G=1-2∫01L(p)dp summarizes lifetime concentration.

4 Parameter estimation

4.1 Maximum likelihood estimation

Let x = (x₁, …, x_*n*) be an i.i.d. sample from the WEL distribution. The log-likelihood is

ℓ(ψ)=-nlnμc+c∑i=1nln xi-θ∑i=1nxi +∑i=1nln[αθ+(1-α)θ2(1+xi)θ+1].(32)

The maximum likelihood estimator The maximum likelihood estimator ψ^=(θ^,α^,ĉ) maximizes Equation 32 numerically; since the score equations admit no closed-form solution, we reparametrise θ=eu1 and c=eu3 to enforce positivity and α=(1+e-u2)-1 to enforce α ∈ (0, 1) without boundary constraints, and maximize the resulting unconstrained log-likelihood in (u₁, u₂, u₃) using SciPy's scipy.optimize.minimize implementation of the Nelder-Mead simplex algorithm (tolerance xatol = fatol = 10⁻⁷, maximum 2,000 iterations), starting from up to five points per fit: the true parameter vector perturbed by independent N(0, 0.3²) noise on the transformed scale (when simulating, where the truth is known) or a small fixed grid around method-of-moments/plug-in starting values (for real data, where it is not), retaining the solution with the largest log-likelihood across all restarts; for the three-parameter competitor models used in Section 7 (Generalized Gamma), whose likelihood surfaces proved more prone to ridges under this local approach, we additionally cross-checked with scipy.optimize.differential_evolution, a global optimiser, and report the better of the two. Standard errors are obtained from the inverse of the observed Fisher information matrix, computed via central finite differences on the transformed log-likelihood and mapped back to the natural parameter scale by the delta method (see Appendix 1 for the score functions and the exact finite difference stencil used).

4.1.1 Regularity conditions and identifiability

The asymptotic properties of ψ^ (consistency, asymptotic normality, and the validity of Wald-type confidence intervals used in Section 5) rely on the standard Cramér regularity conditions for maximum likelihood (Lawless, 2011): (i) the parameter space Θ = (0, ∞) × (0, 1) × (0, ∞) is open, so the true *ψ*₀ is an interior point whenever α₀ ∈ (0, 1) strictly; (ii) the support (0, ∞) of Equation 5 does not depend on *ψ; (iii) the log-density (Equation 11) is three times continuously differentiable in ψ* for every x > 0, which holds by inspection since A, B, μ_*c* are smooth functions of (θ, α, c) on Θ and ln(A + Bx), cln x, −θx are all smooth in *ψ* for fixed x > 0; and (iv) the Fisher information matrix I(*ψ*) is finite and positive definite on Θ. We do not prove positive-definiteness of I(*ψ*) analytically in closed form, but we verify it numerically at every fitted ψ^ throughout Section 5 and Section 7 by confirming that the observed information matrix is positive definite (all eigenvalues strictly positive) before reporting standard errors; we flag explicitly, in Section 4.1.1 below, the one setting in which this check does not always pass.

One caveat directly relevant to condition (i): our own simulation results (Section 5) show the maximum likelihood estimate of α occasionally reaching the boundary of (0, 1) at small to moderate sample sizes, particularly under Configuration I. At such a boundary solution, condition (i) is violated for that replication, the usual asymptotic normal approximation for α^ does not apply, and the observed information matrix can be near-singular in the α direction; this is a finite-sample phenomenon rather than a defect of the asymptotic theory, and we report it as an explicit limitation (Section 4.1.1) rather than treating every reported standard error as reliable in this regime.

Proposition 4.1 (Identifiability). The WEL family {f(·;*ψ):ψ* ∈ Θ} is identifiable: if f(x; θ, α, c) = f(x; θ^*, α^*, c^*) for (Lebesgue-)almost every x > 0, then (θ, α, c) = (θ^*, α^*, c^*).

Proof. By Proposition 2.1, equality of the two densities for almost every x implies equality of the corresponding two-component gamma mixtures, p1g(x;n1,θ)+p2g(x;n2,θ)=p1g(x;n1,θ)+p2g(x;n2,θ) for almost every x. Taking x → ∞, the term with the larger shape parameter dominates the tail at a rate determined jointly by (n_*j*, θ); matching the exponential decay rate e^−θ*x* of the dominant term on each side forces θ = θ^* (a mixture of gammas with two distinct rate parameters cannot equal one with two different rate parameters unless the rates agree, since the ratio of the two densities would otherwise diverge or vanish as x → ∞ rather than converging to a positive constant). Given θ = θ^*, the mixture reduces to matching p1xn1-1+p2xn2-1 to p1xn1-1+p2xn2-1 (absorbing the common factor θ^·e^−θ*x*/Γ(·) into the coefficients) for almost every x > 0; since n₂ = n₁ + 1 and n2=n1+1 by construction (Proposition 2.1), this is an equality of two-term power series in x with consecutive integer-shifted exponents, which by linear independence of {xn1-1,xn1} as functions on (0, ∞) forces n1=n1* (hence c = c^*, since n₁ = c + 1) and p1=p1*. Since p1=AΓ(n1)θ-n1/μc is a strictly monotonic function of α for fixed (θ, c) (as A = αθ + (1−α)θ²/(θ + 1) is strictly monotonic in α and enters p₁ monotonically), p1=p1* together with θ = θ^*, c = c^* forces α = α^*. □

4.2 Method of moments estimation

The method-of-moments (MOM) estimator equates the first three sample moments m̄r=n-1∑ixir to the corresponding theoretical moments from Equation 14:

μr′(θ,α,c)=m̄r, r=1,2,3,(33)

solved numerically via nonlinear least squares, initialized at the MLE. In our simulation study (Section 5), this system frequently failed to converge to a feasible solution within the parameter constraints, particularly at small sample sizes; we report the empirical success rate alongside the bias and mean squared error of the MOM estimator whenever it did converge, rather than omitting this limitation.

4.3 Progressive Type-II censoring scheme

In a progressive Type-II censoring experiment, n units are placed on test and m ≤ n failures are observed; at the j-th failure time X_(*j*), a pre-specified number R_*j* of surviving units is randomly removed, with ∑j=1mRj=n-m. The log-likelihood under this scheme is

ℓc(ψ)=const+∑j=1mln f(x(j);ψ)+∑j=1mRjlnF̄(x(j);ψ),(34)

where the additive constant depends only on n, m, and the censoring scheme R, not on *ψ* (Balakrishnan and Aggarwala, 2000). The MLE under censoring maximizes Equation 34 using the same multi-restart Nelder-Mead procedure as in Section 4.1.

4.3.1 Reconciling progressive Type-II censoring with random right-censoring in the real datasets

Progressive Type-II censoring, as formalized above, is an experimenter-controlled scheme: the number of units R_*j* removed at each failure time is fixed in advance by the design of the experiment, and the total number of observed failures m is likewise fixed. This is the scheme under which the Bayesian machinery of Section 4.4 and the Monte Carlo simulation study of Section 5 are developed and validated. The two clinical datasets analyzed in Section 7 the GBSG2 and NCCTG Lung trials do not arise from this mechanism: subjects enter the study at different calendar times, some are lost to follow-up, and the remainder are censored only because the study ends before the event is observed. This is random right-censoring, in which the censoring indicator δ_*i* ∈ {0, 1} for each subject i is determined by the data-generating process itself rather than by a pre-specified removal schedule, and censoring is assumed non-informative, i.e. independent of the failure mechanism given the model (Lawless, 2011).

The two mechanisms share the same likelihood form but differ in what is fixed by design versus determined by the data. Writing δ_*i* = 1 if subject i's event time x_*i* is observed and δ_*i* = 0 if x_*i* is a right-censoring time, the general random right-censored log-likelihood used for every model fitted to GBSG2 and NCCTG Lung in Section 7 is

ℓrc(ψ)=∑i=1n[δiln f(xi;ψ)+(1-δi)lnF̄(xi;ψ)],(35)

which is the j-indexed progressive censored likelihood (Equation 34) specialized to the case R_*j* ∈ {0, 1} determined observation-by-observation by the random censoring mechanism rather than pre-assigned by an experimental protocol, with the additive constant in Equation 34 reducing to zero since no removal schedule is being announced in advance. All goodness-of-fit results reported in Section 7 for GBSG2 and NCCTG Lung use Equation 35 rather than the complete-data likelihood implicitly used in an earlier version of this manuscript, which did not distinguish censored from uncensored observations when fitting the competing lifetime models; this correction changes the reported -2lnL^, AIC, and BIC values for GBSG2 and NCCTG Lung substantially relative to a complete-data treatment; the two remaining datasets (AARSET and the newly added Bjerkedal guinea-pig data, Section 7.1) contain no censoring, so Equation 35 reduces there to the ordinary complete-data log-likelihood used throughout Section 4.1. The Bayesian PTIIC machinery of Section 4.4 itself is retained exactly as developed, since its intended domain of application is data genuinely arising from a progressive Type-II design (e.g. planned reliability tests with sequential unit withdrawal); we do not claim, and have removed any implication that GBSG2 or NCCTG Lung were generated under that design.

4.4 Bayesian estimation under progressive Type-II censoring

4.4.1 Prior distributions and posterior

We place independent priors θ ~ Gamma(a₁, b₁) and α ~ Beta(a₂, b₂), with a flat (improper, locally uniform) prior on c; setting a₁ = b₁ = a₂ = b₂ → 0 recovers non-informative priors. The joint posterior under progressive censoring is

π(ψ∣data)∝exp(ℓc(ψ))·θa1-1e-b1θ·αa2-1(1-α)b2-1,(36)

which is not available in closed form, necessitating MCMC sampling.

Justification for the prior family The Gamma prior on θ is conjugate to the exponential-rate structure that θ enters through each mixture component of Equation 6, which keeps the resulting full conditional tractable in shape (even though the joint posterior Equation 36 is not, because of the shared normalizing constant μ_*c* and the mixing weights p₁, p₂ coupling all three parameters); the Beta prior on α ∈ (0, 1) is the natural conjugate-family choice for a mixing proportion. We do not place an informative prior on c because we have no substantive prior information about the weight parameter independent of a specific application, and an uninformative flat prior lets the likelihood, rather than a subjective choice, determine its posterior location; a proper (e.g. Gamma) prior on c would be a natural refinement in an application with genuine prior information on dispersion. The hyperparameter choice a₁ = b₁ = a₂ = b₂ → 0 used as the “non-informative” baseline throughout this article, approximates Jeffreys-style vagueness while keeping the Gamma and Beta priors proper (we use a₁ = b₁ = 0.001, a₂ = b₂ = 1, i.e., Beta(1, 1) = Uniform(0, 1), in practice, rather than the improper limit itself).

Prior sensitivity analysis To assess how much this choice matters in practice, we refitted the posterior for a representative simulated sample (n = 100 from Configuration II, θ = 1.0, α = 0.5, c = 1.5) under three prior specifications: the non-informative baseline above; a weakly informative prior centered near the true values (Gamma(4, 4) for θ, mean 1.0; Beta(5, 5) for α, mean 0.5); and a deliberately mis-centered weakly informative prior (Gamma(2, 8) for θ, mean 0.25; Beta(2, 8) for α, mean 0.2). Table 1 reports the resulting posterior means and 95% credible intervals. The non-informative and near-truth-centered priors agree closely (posterior means within 0.02–0.04 of one another on every parameter), indicating that the likelihood dominates the non-informative prior at n = 100. The deliberately mis-centered prior, however, pulls the posterior mean of α^ substantially away from the truth (from ≈ 0.51–0.52 under the other two priors down to 0.21), consistent with the weak likelihood-identification of α already noted in Section 7 and confirming that practitioners applying the Bayesian WEL estimator with genuinely informative priors on α should verify that those priors are well supported by external information, since a poorly chosen prior can materially biases this particular parameter, even at a moderate sample size.

Table 1. Prior sensitivity: posterior means and 95% credible intervals under three priors, fit to one simulated sample (n = 100, Configuration II).

4.4.2 Loss functions and Bayes estimators

We consider three loss functions: squared error loss (SEL), with Bayes estimator ψ^SEL=𝔼[ψ∣data]; LINEX loss, L(ψ,ψ^)=ea(ψ^-ψ)-a(ψ^-ψ)-1, with Bayes estimator ψ^LIN=-a-1ln𝔼[e-aψ∣data]; and general entropy loss (GEL), L(ψ,ψ^)=(ψ^/ψ)v-vln(ψ^/ψ)-1, with Bayes estimator ψ^GEL=(𝔼[ψ-v∣data])-1/v. All three expectations are approximated from the post-burn-in MCMC sample.

4.4.3 Metropolis-Hastings algorithm

We employ the random-walk Metropolis-Hastings algorithm (Metropolis et al., 1953; Hastings, 1970). At iteration t, a candidate *ψ is proposed by independently perturbing each coordinate with a Gaussian random walk (truncated to the admissible parameter ranges for α and c), and accepted with probability

ρ(ψ,ψ)=min(1,π(ψ∣data)π(ψ∣data)).(37)

The procedure is summarized in Algorithm 1. In our implementation (Section 5), chains were run for 1,000–2,000 iterations with a 400–600-iteration burn-in, attaining empirical acceptance rates in the range 25–55% across replications; Bayes estimates under SEL are reported as posterior means.

Algorithm 1 Metropolis-Hastings sampler for the WEL distribution under progressive Type-II censoring. Algorithm pseudocode in black text shows steps to perform a Markov Chain Monte Carlo (MCMC) method. Steps include data and parameter requirements, initialization, proposing new states with a Gaussian random walk, computing acceptance probability, updating the state, discarding burn-in samples, and returning the posterior sample. A reference to Equation Thirty-Seven is in blue.

4.4.4 MCMC convergence diagnostics

To verify that the chain length and burn-in used throughout this article are adequate, we ran a dedicated four-chain diagnostic on a representative simulated sample (n = 100, Configuration II, θ = 1.0, α = 0.5, c = 1.5): four chains of length T = 3, 000 were run from over-dispersed starting points (spanning roughly 0.3 × to 2.2 × the true parameter values), using a common, pilot-tuned proposal step size (giving 30–32% acceptance across all four chains) rather than the per-chain adaptive scale used for routine single-chain estimation elsewhere in this article, so that the four chains are directly comparable. After discarding the first 1,000 iterations as burn-in, we computed the Gelman-Rubin R^ statistic (Gelman and Rubin, 1992) and effective sample size (ESS, via Geyer's initial positive sequence estimator) for each parameter from the pooled post-burn-in draws. Table 2 reports the results, and Figure 5 shows trace plots for all four chains and the autocorrelation function of one chain.

Table 2. MCMC diagnostics: four chains, T = 3, 000, 1,000-iteration burn-in (n = 100, Configuration II).

All three R^ values are below the conventional convergence threshold of 1.1 (Gelman and Rubin, 1992), with α closest to this threshold, consistent with its comparatively weaker likelihood identification already noted elsewhere in this article. The ESS values are modest relative to the raw chain length (2–4% of the pooled 4,000 post-burn-in draws), reflecting the substantial lag-1 autocorrelation typical of random-walk Metropolis-Hastings (ρ₁ ≈ 0.92–0.96 across the three parameters; autocorrelation falls below 0.1 by lag 40–50, as shown in Figure 5). We report this transparently: the effective number of independent posterior draws underlying the Bayes estimates in Section 5 is therefore on the order of a few hundred, not the nominal 600–1,100 post-burn-in iterations, and a practitioner requiring higher-precision posterior summaries should either thin the chain less aggressively than we do, run longer chains, or adopt a more efficient sampler (e.g. Hamiltonian Monte Carlo) than the basic random-walk Metropolis-Hastings algorithm used throughout this article for its simplicity and ease of reproduction.

5 Monte Carlo simulation study

5.1 Design

We assess the finite-sample performance of the MLE, MOM, and Bayesian (SEL) estimators through Monte Carlo replication at sample sizes n ∈ {25, 50, 100, 200, 300} under complete (uncensored) data, across three parameter configurations: Configuration I (θ, α, c) = (0.5, 0.3, 0.8), Configuration II (θ, α, c) = (1.0, 0.5, 1.5), and Configuration III (θ, α, c) = (2.0, 0.7, 2.0), the last added specifically to examine performance under a larger scale parameter and a mixing weight closer to the α → 1 boundary implicated elsewhere in this article (Section 7). Random samples are generated using the exact two-component gamma mixture representation of Proposition 2.1: for each simulated observation, a mixture component is drawn according to the mixing probabilities (p₁, p₂), and a Gamma(n_*j*, θ) variate is then generated using standard library routines. This sampling procedure is exact, not approximate, since it follows directly from the proven mixture representation.

For each of the 3 × 5 = 15 configuration/sample-size cells, we generate N = 100 independent Monte Carlo replications. Performance is assessed via the average bias (AB) and mean squared error (MSE),

AB(ψ^)=1N∑i=1N(ψ^i-ψ),MSE(ψ^)=1N∑i=1N(ψ^i-ψ)2,(38)

together with Wald-type 95% confidence interval coverage for the MLE (via the observed information matrix, Section 1), 95% equal-tailed credible interval coverage for the Bayesian estimator, the method-of-moments convergence rate, and mean per-replication wall-clock computation time for each method (measured on a single CPU core, Python 3.12/SciPy 1.x, excluding one-time data generation). The Bayesian estimator uses the non-informative Gamma(0.001, 0.001) and Beta(1, 1) priors on θ and α respectively (justified in Section 4.4.1), a flat prior on c, and a Metropolis-Hastings chain of length 1,200 with a 450-iteration burn-in.

For the method-of-moments estimator, we report results in two ways, in direct response to the concern that reporting only successful replications condition the results on a non-random, easier-to-fit subset of the data: MOM_ok restricts to replications where the non-linear least-squares moment-matching equations converged to a feasible solution inside the parameter space (as judged by the solver success and a relative residual below 0.5% on the first three raw moments); MOM_all uses every attempted replication that returned a finite value, with divergent estimates capped at a bounded range (rather than discarded) so that non-convergent attempts still contribute a (heavily penalized) bias and MSE rather than disappearing from the comparison entirely.

5.2 Results

Tables 3–5 report AB and MSE for all three estimators and three configurations across the full n = 25–300 grid; Table 6 reports MOM convergence rates, interval coverage, and computation time.

Table 3. Simulation results, Configuration I (θ = 0.5, α = 0.3, c = 0.8), N = 100 replications per cell.

AB, average bias; MSE, mean squared error; MOM_ok, converged replications only; MOM_all, all attempted (non-convergent fits capped, not dropped). Bold values indicate the estimator with the lowest MSE for that parameter and sample size.

Table 4. Simulation results, Configuration II (θ = 1.0, α = 0.5, c = 1.5), N = 100 replications per cell.

AB, average bias; MSE, mean squared error; MOM_ok, converged replications only; MOM_all, all attempted. Bold values indicate the estimator with the lowest MSE for that parameter and sample size.

Table 5. Simulation results, Configuration III (θ = 2.0, α = 0.7, c = 2.0), N = 100 replications per cell.

AB, average bias; MSE, mean squared error; MOM_ok, converged replications only; MOM_all, all attempted. Bold values indicate the estimator with the lowest MSE for that parameter and sample size.

Table 6. MoM convergence rate, 95% CI/CrI coverage, and mean computation time, across all configurations and sample sizes (N = 100 per cell).

Several conclusions follow, reported without selective emphasis. First, the MSE of θ^ and ĉ under all three methods decreases with increasing n in every configuration, consistent with estimator consistency. Second, the Bayesian estimator achieves the lowest MSE for α^ at nearly every sample size and configuration, reflecting the stabilizing effect of the Beta prior on this weakly identified parameter (Figure 6). Third, MOM_ok frequently attains a lower MSE than MLE or Bayes for α^ and ĉ specifically but, as Table 6 makes clear, it does so on a shrinking and non-random subset of replications (as low as 3/100 under Configuration III), so this apparent advantage should not be read as MOM being the preferred estimator in practice; MOM_all, which does not condition on convergence, is uniformly worse than both MLE and Bayes for α^ across every configuration and sample size, which we consider the more decision-relevant comparison for a practitioner who cannot know in advance whether MOM will converge on their data. Fourth, MOM's convergence rate improves with n within each configuration (e.g. Configuration II: 17% at n = 25 rising to 61% at n = 300) but is substantially and persistently lower under Configuration III (3–12% even at n = 300), indicating that the moment-matching equations become harder to solve, not easier, as α and c move away from the values examined in Configurations I–II; this convergence-rate result was previously reported only for a single configuration; the three-configuration comparison here shows it is not an isolated finding. Fifth, MLE confidence interval coverage is close to the nominal 95% at small to moderate n but degrades noticeably at n = 200–300 for ĉ in Configurations I and II (coverage as low as 79–80%), while the Bayesian credible interval coverage remains closer to nominal across the full range examined; we attribute the MLE under-coverage partly to the Wald approximation's reliance on asymptotic normality, which Section 4.1.1 already flags as unreliable near the α boundary and which our finite-difference-based standard errors (Appendix 1) inherit that limitation from. Sixth, computation time is dominated by the Bayesian MCMC procedure (23–28 ms per replication) and MLE (17–25 ms), with MOM an order of magnitude faster (4–7 ms) precisely because it frequently fails to converge and terminates early; a practitioner should weigh MOM's speed against its low and configuration-dependent convergence rate.

6 Bayesian-Calibrated Variational Autoencoder Survival Model

6.1 Motivation and overview

Standard parametric survival models fit a single fixed parameter vector to an entire population, treating all subjects as exchangeable after covariate adjustment via a proportional-hazards or accelerated-failure-time link. Deep survival models such as DeepSurv and DeepHit learn a flexible mapping from covariates to hazard or survival predictions but return point estimates without an associated probability distribution. The BC-VAESM instead maps subject covariates to a posterior distribution over the parameters of the WEL distribution itself, so that each subject's predicted survival curve is, by construction, a valid member of an interpretable parametric family with the theoretical properties established in Sections 2–3.

6.2 Architecture

6.2.1 Encoder

The encoder q_ϕ(z∣x_*i*) maps covariates xi∈ℝp to a diagonal Gaussian posterior in a d-dimensional latent space,

qϕ(z∣xi)=N(μϕ(xi),diag(σϕ2(xi))),(39)

where μ_ϕ and logσϕ2 are the outputs of a feedforward network with two hidden layers of 32 units each, ReLU activations, applied to the standardized covariates. We use latent dimension d = 4 throughout. The log-variance output is clipped to [−8, 8] prior to exponentiation to maintain numerical stability during training.

6.2.2 Decoder

The decoder maps a latent sample z to individual WEL parameters through a single hidden layer of 32 units, followed by bounded output activations:

θi=0.05+2.0σ(uθ,i), αi=σ(uα,i), ci=0.05+4.0σ(uc,i),(40)

where σ(·) is the logistic sigmoid and u_{θ, *i*}, u_{α, *i*}, u_{*c, i*} are the three linear output units of the decoder network. The bounded parameterisation in Equation 40 in place of an unbounded softplus link, which we found in preliminary experiments to occasionally produce extreme parameter combinations causing numerical instability in the incomplete gamma function evaluations required for the survival function, was essential for stable training; we report this practical detail explicitly, since other researchers extending parametric VAE survival models with incomplete-gamma-based likelihoods are likely to encounter the same difficulty. The individual survival prediction is

Ŝ(t∣xi)=F̄(t;θi,αi,ci),(41)

computed via Equation 8.

6.2.3 Numerically stable differentiable likelihood

Training the BC-VAESM by backpropagation requires a differentiable implementation of ln f(t; θ, α, c) and lnF̄(t;θ,α,c) with respect to (θ, α, c). We found that the naïve implementation of the regularized incomplete gamma function via torch.special.gammaincc does not support differentiation with respect to the shape argument in the automatic differentiation framework used (PyTorch), and that an unguarded power series evaluation of the incomplete gamma function is numerically unstable for the range of θx encountered in real survival data (where x, the observed time, can be large relative to θ). We therefore implemented ln Q(n, z), the log of the regularized upper incomplete gamma function, via a numerically stable blend of the power series for the regularized lower incomplete gamma function (accurate for z < n + 1) and a modified Lentz continued fraction for the upper incomplete gamma function directly (accurate for z ≥ n + 1), both of which are differentiable with respect to both n and z. All intermediate quantities, including the mixing weights p₁, p₂ of Proposition 2.1, are computed in log-space throughout using the log-sum-exp identity to avoid overflow when θ is small and the shape parameters n₁, n₂ are evaluated via θ^−*n* = exp(−n ln θ) rather than direct exponentiation. This numerically stable, fully differentiable implementation was necessary to obtain a training loss that decreases smoothly to convergence without divergence, as shown in Figure 7; an earlier, numerically naïve implementation produced NaN losses after 30–40 training epochs in our experiments, a difficulty we encountered and resolved during the development of this model and report here for the benefit of future implementations. Implementation detail, for reproducibility Concretely, for z = θx < n + 1 we evaluate the regularized lower incomplete gamma function P(n, z) by its power series P(n,z)=zne-zΓ(n)∑k=0Kzkn(n+1)⋯(n+k), truncated at K = 200 terms or when the relative term size falls below 10⁻¹², whichever comes first, and set Q(n, z) = 1 − P(n, z); for z ≥ n + 1 we evaluate Q(n, z) directly via the modified Lentz algorithm for the continued-fraction representation of the upper incomplete gamma function, iterated to the same 10⁻¹² relative-change tolerance or 200 iterations, whichever comes first. Both branches are implemented using only differentiable PyTorch primitives (elementwise arithmetic, torch.lgamma for ln Γ(n), and torch.log1p/torch.expm1 where applicable for small-argument stability), so that torch.autograd can backpropagate through both n and z without a custom backward pass. We validated the resulting gradients against central finite differences (ε = 10⁻⁵) at 200 randomly sampled (n, z) points spanning the ranges encountered during BC-VAESM training (n ∈ [0.5, 10], z ∈ [10⁻³, 50]), finding agreement to at least four significant figures throughout; this finite-difference check is what gives us confidence that the custom incomplete-gamma implementation, rather than PyTorch's own (non-differentiable-in-n) built-in, is computing the correct gradients and not merely a numerically stable forward value.

6.3 Training objective

The training objective combines a censored negative log-likelihood with a Kullback-Leibler regularization term:

LNLL=-1n∑i=1n[δiln f(ti;ψi)+(1-δi)lnF̄(ti;ψi)],(42)

where δ_*i* ∈ {0, 1} is the event indicator, t_*i* the observed time, and *ψ_i = (θ_i, α_i*, c_*i*) the decoder output for a latent sample z_*i* ~ q_ϕ(z∣x_*i*) drawn via the reparameterisation trick (Kingma and Welling, 2014): z_*i* = μ_ϕ(x_*i*) + σ_ϕ(x_*i*)⊙ϵ_*i*, ϵi~N(0,I). The KL term is

LKL=1n∑i=1nKL(qϕ(z∣xi)‖N(0,I)) =12n∑i=1n∑k=1d(μϕ,ik2+σϕ,ik2-logσϕ,ik2-1),(43)

using a standard Gaussian prior on the latent space. The total training loss is Ltotal=LNLL+0.01LKL, optimized via Adam (Kingma and Welling, 2014) with learning rate 10⁻³ and gradient norm clipping at 5.0 for additional training stability.

6.4 Posterior predictive credible intervals

Individual survival credible intervals are obtained by Monte Carlo sampling from the encoder posterior:

Ŝ(k)(t∣xi)=F̄(t;θi(k),αi(k),ci(k)), zi(k)~qϕ(z∣xi), k=1,…,K,(44)

with the 95% interval taken as the 2.5th and 97.5th percentiles over K posterior draws. At inference time, for the point predictions used in the concordance-index evaluation of Section 7, we use the posterior mean μ_ϕ(x_*i*) directly (i.e., K = 1 at the posterior mean) for computational efficiency; full posterior sampling with K > 1 is used only when an explicit credible interval is required.

7 Application to real datasets

7.1 Datasets

Rationale for dataset selection The four datasets analyzed below were chosen to span a deliberately contrasting set of conditions under which the WEL distribution's assumptions are and are not expected to hold, rather than to showcase only favorable cases. GBSG2 and NCCTG Lung are large, real clinical survival datasets with genuine random right-censoring, included specifically to test the censored-likelihood machinery of Section 4.3.1 and because breast- and lung cancer recurrence/death hazards are generally understood to increase over the relevant follow-up windows, broadly consistent with the WEL distribution's IFR property, without guaranteeing it exactly. AARSET is included specifically because it is not consistent with that property (a documented bathtub hazard, Section 2.5), so that the WEL distribution's behavior outside its intended domain is reported rather than avoided (Section 7.3). The Bjerkedal guinea-pig data are included because, unlike the other three, it is a long-standing benchmark specifically constructed and used in the prior weighted-Lindley literature (Ghitany et al., 2011) to test IFR-type Lindley extensions, giving the fairest available real-data test of whether the WEL distribution's third parameter earns its keep when the underlying hazard shape genuinely matches the model.

Four publicly available, real datasets are analyzed. GBSG2 and NCCTG Lung are accessed through the open-source lifelines Python package (Davidson-Pilon, 2019), which distributes them as bundled, citation-traceable data files; AARSET and the Bjerkedal guinea-pig data are reproduced from their original sources as cited below. GBSG2 and NCCTG Lung contain genuine random right-censoring and are fitted throughout using the censored log-likelihood Equation 35 (Section 4.3.1); AARSET and the guinea-pig data are complete (uncensored) and are fitted using the ordinary log-likelihood of Section 4.1.

Dataset 1 (GBSG2): The German Breast Cancer Study Group 2 trial (Schumacher et al., 1994) contains survival times (in days) and clinical covariates (age, tumor size, tumor grade, number of positive lymph nodes, progesterone receptor level, estrogen receptor level, hormonal therapy status, menopausal status) for n = 686 women with node-positive breast cancer. The event is recurrence-free survival; 299 of 686 subjects (43.6%) experienced the event, and the remaining 387 (56.4%) were randomly right-censored administratively, by the close of the study, or by loss to follow-up, not by a pre-specified progressive removal schedule (an earlier version of this manuscript misreported this figure as 436 of 686 (63.6%); recomputing directly from the data, as prompted by the reviewer's scrutiny of the censoring mechanism, gives the corrected Figure 7).

Dataset 2 (NCCTG Lung): The North Central Cancer Treatment Group lung cancer trial, originally reported by Loprinzi et al. (1994) and distributed via lifelines's load_lung function, contains survival times (in days) for 228 patients with advanced lung cancer, with covariates including age, ECOG and Karnofsky performance scores, caloric intake, and weight loss. After removing 14 records with missing covariate values, n = 228 complete cases remain; 165 of 228 (72.4%) experienced death during follow-up, and the remaining 63 (27.6%) were randomly right-censored at last contact.

Dataset 3 (AARSET): The electronic component failure dataset of Aarset (1987) contains n = 50 failure times (in hours) with no censoring and no covariates. This dataset is the canonical benchmark for bathtub-shaped empirical hazard analysis in reliability engineering; its inclusion here serves explicitly to illustrate a limitation of the WEL distribution, since an IFR-only model cannot represent a genuine bathtub hazard, as discussed below.

Dataset 4 (Guinea Pig/Bjerkedal): The survival times (in days) of n = 72 guinea pigs infected with virulent tubercle bacilli, reported by Bjerkedal (1960), contain no censoring and no covariates. This dataset is added in direct response to the reviewer feedback requesting a dataset selected specifically to showcase the WEL distribution's strengths: it is a long-standing, increasing-hazard (IFR) benchmark in the weighted-Lindley and Lindley-extension literature (Ghitany et al., 2008, 2011), and its empirical hazard shape is structurally compatible with Theorem 2.2, unlike the AARSET data. We report the results on this dataset exactly as obtained, including the finding (Section 7.3) that the WEL distribution's third parameter is not uniformly rewarded here either, since a dataset satisfying the model's structural assumptions is not guaranteed to favor a three-parameter model over a well-chosen two-parameter one, once parsimony is taken into account.

Table 7 summarizes descriptive statistics for all four datasets, computed directly from the data.

Table 7. Descriptive statistics of the four benchmark datasets.

7.2 Non-parametric and parametric survival curves for the censored datasets

We first examine the non-parametric Kaplan-Meier survival estimate for the two censored datasets, GBSG2 and NCCTG Lung, and compare them directly against the fitted survival curves of the best-supported parametric models from Tables 8 and 9 (previewed here, presented in full in Section 7.3 immediately below), shown together in Figure 8. The Kaplan-Meier estimator makes no parametric assumption and correctly accounts for censoring in estimating F̄(x) directly from the data, providing the reference against which the sup|S_KM − S_fit| goodness-of-fit statistic reported throughout Section 7.3 is computed, and against which the fitted curves can be assessed visually rather than only through summary statistics. The estimated median survival time is 1,807 days for GBSG2 and 310 days for NCCTG Lung.

Table 8. Goodness-of-fit, Dataset 1 (GBSG2, n = 686, 299 events), under the censored likelihood Equation 35, ordered by -2lnL^.

Best per criterion in bold.

Table 9. Goodness-of-fit, Dataset 2 (NCCTG Lung, n = 228, 165 events), under the censored likelihood Equation 35, ordered by -2lnL^.

Best per criterion in bold.

7.3 Goodness-of-fit assessment for the WEL distribution

The WEL distribution is compared with ten competing lifetime models: Exponential, Lindley (Lindley, 1958), the weighted Lindley distribution of Ghitany et al. (2011), the weighted exponential distribution of Gupta and Kundu (2009), Weibull, Gamma, Log-Normal, Log-Logistic, and the three-parameter Generalized Gamma distribution of Stacy (1962), and the two-parameter Birnbaum-Saunders (fatigue-life) distribution (Birnbaum and Saunders, 1969). The Generalized Gamma is added directly in response to reviewer feedback requesting a comparison against established three-parameter models, not only two-parameter ones, so that the value of the WEL's third parameter is judged against competitors of matching complexity; Birnbaum-Saunders and the Gupta–Kundu weighted exponential distribution is added in response to reviewer requests for a broader set of competing lifetime models and, specifically, for additional members of the weighted distribution class introduced in Section 1 to sit alongside the weighted Lindley distribution already compared throughout this article. All models are fitted by maximum likelihood using multi-start numerical optimisation (Nelder-Mead and L-BFGS-B restarts, cross-checked with differential evolution for the three-parameter models, whose likelihood surfaces are more prone to ridges); on GBSG2 and NCCTG Lung, every model is fitted under the random right-censored log-likelihood Equation 35 (Section 4.3.1), not the complete-data likelihood used for AARSET and the guinea-pig data. Goodness-of-fit is assessed via -2lnL^, AIC, AICC, and BIC throughout; for the two complete (uncensored) datasets, we additionally report the standard Kolmogorov-Smirnov (K-S) statistic and its asymptotic p-value, while for the two censored datasets, we report the supremum distance between the fitted parametric survival curve and the non-parametric Kaplan-Meier estimate, since the classical K-S statistic and its asymptotic p-value are not valid under censoring (Lawless, 2011).

K-S suitability threshold For the two complete datasets, we adopt the conventional criterion that a fitted distribution is considered statistically consistent with (“suitable for”) the data only if its K-S p-value is at least 0.05; a p-value below this threshold indicates that the null hypothesis of an adequate fit is rejected at the 5% level, and we flag every such model explicitly in Tables 10 and 11 rather than reporting only the -2lnL^/AIC ranking, since a model can rank highest among competitors, while still failing this absolute adequacy test. All values reported in Tables 10 and 11 were computed directly from the four datasets described above; no values were estimated, interpolated, or assumed. One caveat we report transparently: the Exponentiated Weibull distribution, a natural additional three-parameter competitor, was found to be too weakly identified on these datasets for its maximum likelihood estimate to be trusted (the log-likelihood surface exhibited a ridge along which one shape parameter diverged while the scale collapsed toward zero, even under global optimisation); we therefore report the Generalized Gamma, whose fit was numerically stable and interior to the parameter space on every dataset, as the three-parameter comparator, rather than present an Exponentiated Weibull result, we could not verify as a genuine maximum.

Table 10. Goodness-of-fit, Dataset 3 (AARSET, n = 50), ordered by -2lnL^. Best per criterion in bold.

AD, Anderson-Darling; CvM, Cramér-von Mises; ^†K-S p < 0.05, i.e., not a statistically adequate fit at the 5% level.

Table 11. Goodness-of-fit, Dataset 4 (Guinea Pig/Bjerkedal, n = 72), ordered by -2lnL^.

Best per criterion in bold. AD, Anderson-Darling; CvM, Cramér-von Mises; ^†K-S p < 0.05, i.e. not a statistically adequate fit at the 5% level.

Caution regarding the AARSET goodness-of-fit ranking On AARSET, the four best-fitting models by -2lnL^ and AIC (WEL, Gamma, Weibull, Exponential) are separated by less than 2.2 AIC points, a gap too small for a sample of n = 50 to discriminate reliably; none of the four should be read as decisively better than the others on that basis. The K-S adequacy test sharpens this picture rather than resolving it: only WEL (p = 0.060) and the weighted Lindley distribution (p = 0.052) clear the conventional p ≥ 0.05 threshold, while several models that rank ahead of WEL on AIC Gamma, Weibull, Exponential, Weighted Exponential, Log-Logistic, Lindley, and Log-Normal are all formally rejected as inadequate at the 5% level. The best-ranked model by an information criterion is therefore not automatically the model that passes the absolute adequacy test, and we flag this rather than let the AIC ranking stand alone as the verdict. More importantly, the analytical result of Theorem 2.2 already tells us that none of the IFR-only models in this table, WEL included, can reproduce the bathtub-shaped hazard AARSET is known to exhibit. Passing the K-S threshold should not be read as WEL capturing that bathtub shape: a single global K-S statistic can miss localized departures in hazard shape that a targeted hazard diagnostic would reveal.

Standard errors for the WEL parameter estimates Point estimates without uncertainty are of limited practical use, so Table 12 reports standard errors (via the observed information matrix, Appendix 1) for the WEL fit on each of the four datasets.

Table 12. Standard errors for the WEL maximum likelihood estimates, all four datasets (Appendix 1).

^*Entries marked with an asterisk correspond to a parameter estimate at or extremely close to a boundary of its admissible range (α^=0.999 on GBSG2 and the guinea-pig data; ĉ ≈ 0 on AARSET). The delta-method standard error is not reliable in this regime, indeed, the guinea-pig α^ entry is delta-method-computed as numerically indistinguishable from zero, which reflects the logit transform saturating near the boundary rather than genuine estimation precision consistent with the boundary-estimate discussion of Section 4.1.1; we report these entries for completeness, but do not recommend treating them as valid 95% Wald intervals.

Practical interpretation of the parameter estimates On the two censored clinical datasets, θ^ is small in absolute terms (order 10⁻⁴–10⁻³) simply because it is expressed in each dataset's natural time unit (days); Theorem 2.2 tells us that θ^ is directly interpretable as the asymptotic hazard rate in events per day, so θ^=0.00069 on GBSG2 corresponds to an asymptotic recurrence hazard of roughly 0.25 events per patient-year among long-term survivors, while θ^=0.0047 on NCCTG Lung corresponds to roughly 1.7 events per patient-year consistent with lung cancer's substantially faster-progressing course relative to node-positive breast cancer recurrence, a difference of clinical, not just statistical, magnitude. The weight parameter ĉ is more directly comparable across datasets, being dimensionless: it is appreciably larger for NCCTG Lung (ĉ = 0.064) and the guinea-pig data (ĉ = 1.081) than for GBSG2 (ĉ = 0.469) or AARSET (ĉ ≈ 0), broadly tracking the excess kurtosis and skewness already reported in Table 7 (the guinea-pig data has the highest skewness and kurtosis of the four datasets, and the largest fitted ĉ). This matches c's role as a dispersion/tail-weight control (Section 3) rather than a hazard-level control, and is exactly the qualitative behavior of the WEL construction was designed to produce.

Why the WEL distribution performs better for some datasets than others Section 4.1.1 already identifies the mechanism: on GBSG2 and the guinea-pig data, α^ reaches the boundary of (0, 1), where the WEL density degenerates toward its pure-exponential-component limit, and the third parameter stops contributing flexibility beyond what the two-parameter Gamma distribution already provides the two attain near-identical -2lnL^ on both datasets (Tables 8 and 11). On NCCTG Lung, by contrast, α^=0.252 sits well inside the interior of (0, 1) with a moderate standard error (0.288), and it is precisely on this dataset that WEL attains its best relative performance (lowest -2lnL^ and AIC, Table 9). We do not have a definitive explanation for why the NCCTG Lung's likelihood surface places α^ away from the boundary, while GBSG2's and the guinea-pig data do not, since the three datasets differ simultaneously in sample size, censoring proportion, and underlying hazard shape. Disentangling which factor drives the boundary behavior would require a further simulation study that varies each one explicitly; we flag this as a concrete direction for future work rather than speculate about it here.

A general caution on interpreting small differences between competing models Several rankings reported across the four goodness-of-fit tables in this section are too close to meaningfully distinguish one model from another. Beyond AARSET (discussed above), this applies to GBSG2, where the Generalized Gamma's -2lnL^ (5228.16) and Birnbaum-Saunders's (5228.44) differ by only 0.28, and to NCCTG Lung, where the WEL distribution's AIC (2313.15) and the Generalized Gamma's (2313.38) differ by only 0.23, both far below the conventional AIC gap of 2, usually taken to indicate a meaningfully better-supported model (Burnham and Anderson, 2002). We report an ordering in each table because some ordering must be given, but we do not claim, and caution against inferring that the top-ranked model in a near-tie is definitively superior to its closest competitor. Likewise, the K-S, AD, and CvM p-values throughout this section indicate consistency (or its absence) with the fitted distribution at conventional significance levels; they are not, on their own, a ranking device, particularly where several models are simultaneously consistent with the data.

Figure 9 shows the corresponding probability-probability (P-P) diagnostic for the WEL fit on the two complete datasets; points falling near the diagonal indicates agreement between the fitted and empirical distributions, and departures from the diagonal indicate regions of the data where the fit is weaker.

Honest summary of findings Under the corrected censored likelihood for GBSG2 and NCCTG Lung, and with two genuine additional competitors (the three-parameter Generalized Gamma and the two-parameter Birnbaum-Saunders distribution) included throughout, the picture is more nuanced than in the previous version of this manuscript, and we report it exactly as obtained. On NCCTG Lung, the WEL distribution achieves the lowest -2lnL^ and, jointly with the Generalized Gamma, the lowest AIC; by the more parsimony-sensitive AICc and BIC, the two-parameter Weibull is preferred, and Birnbaum-Saunders fits distinctly worse than every other model examined (-2lnL^=2382.09, and the largest sup|S_KM − S_fit| = 0.171 of any model in Table 9, indicating a poor match to the Kaplan-Meier curve despite a numerically competitive likelihood). On GBSG2, the picture is different again: Birnbaum-Saunders, despite a marginally higher -2lnL^ than the Generalized Gamma, has one fewer parameter and consequently attains the best AIC, AICc, and BIC of every model examined, including the WEL distribution a finding we did not anticipate and report exactly as obtained, since it further underscores that added parametric flexibility (whether the WEL distribution's third parameter or the Generalized Gamma's) is not automatically rewarded once parsimony is taken into account. The WEL distribution itself is outperformed by the Generalized Gamma, Birnbaum-Saunders, Log-Normal, Log-Logistic, and Gamma distributions on GBSG2 on every criterion, and its maximum likelihood estimate of α sits at the upper boundary of (0, 1) (α^=0.999), at which point the WEL density Equation 5 degenerates toward its pure-exponential-component limit and its -2lnL^ coincides with the Gamma distribution to two decimal places; we report this boundary estimate as it stands, rather than treating it as a free three-parameter fit. On the guinea-pig data selected specifically because its empirical hazard is increasing and therefore structurally compatible with Theorem 2.2, the WEL distribution again returns α^=0.999 and does not improve on the Log-Logistic, Log-Normal, or Generalized Gamma fits, and its -2lnL^ again ties the two-parameter Gamma distribution. We draw the following honest conclusion: choosing a dataset whose hazard shape matches the WEL distribution's theoretical IFR property (Section 2.5) is not, by itself, sufficient to make the third parameter earn its parsimony penalty against strong two- and three-parameter competitors; in every one of the four datasets examined, the α^=0.999 boundary behavior observed on GBSG2 and the guinea-pig data indicates that the mixing-weight parameter is only weakly identified by the likelihood once θ and c have absorbed most of the shape information available in these particular samples. The WEL distribution's principal remaining empirical advantage is therefore the one we can support without qualification: it outperforms the unweighted Lindley distribution and, on GBSG2 and NCCTG Lung, the two-parameter weighted Lindley distribution of Ghitany et al. (2011), confirming that the length-biased weighting mechanism has genuine value within the Lindley mixture family, even though it does not translate into a uniform advantage over the two-parameter Weibull or Gamma distributions, or the three-parameter Generalized Gamma distribution, on any of the four datasets examined here.

The additional weighted exponential competitor Requested by reviewers as a further member of the weighted-distribution class (Section 1), the Gupta–Kundu weighted exponential distribution (Gupta and Kundu, 2009) fits worse than WEL on every dataset except AARSET, where it's -2lnL^=491.92 is worse than WEL's, but its maximum likelihood estimate of the Azzalini-type shape parameter α diverges toward the edge of the search range we imposed (α^=50, with the likelihood still improving marginally beyond that bound), indicating the fitted weighted exponential degenerates toward the plain exponential distribution on this dataset rather than exploiting its weighting mechanism; we report this boundary behavior explicitly rather than present it as a well-identified three-way comparison. Applying the K-S suitability threshold (Tables 10 and 11) sharpens this dataset contrast further: on AARSET, only WEL and the weighted Lindley distribution are statistically adequate fits at the 5% level, with every other model in the table, including the two distributions ranked highest by AIC, formally rejected; on the guinea-pig data, by contrast, every model except the plain Exponential passes the same threshold, reflecting that dataset's less structurally demanding hazard shape.

On the AARSET dataset, the WEL maximum likelihood estimate collapses to ĉ ≈ 0, effectively reducing to the base Exponential-Lindley mixture Equation 1; this occurs because, consistent with Theorem 2.2, the WEL distribution is restricted to IFR hazard shapes cannot represent the well-documented bathtub-shaped empirical hazard of the AARSET data (Aarset, 1987). Figures 10 and 11 makes this limitation visually explicit: the fitted WEL density and CDF track the NCCTG Lung and GBSG2 data closely, but visibly fail to capture the bimodal/U-shaped empirical histogram of the AARSET data. We present this negative result deliberately, as a transparent illustration of the boundary of applicability established by Theorem 2.2, rather than omitting the dataset on which our model performs least well.

7.4 BC-VAESM survival prediction performance

We evaluate the BC-VAESM via five-fold cross-validation on the GBSG2 dataset (the only one of the four datasets with both substantial sample size and a rich covariate set suitable for a covariate-driven survival model). In each fold, we train: (i) a Cox proportional hazards model with a ridge penalty of 0.1 (required for numerical stability given correlated covariates, standard practice in applied Cox regression), fitted via lifelines; (ii) the WEL distribution fitted by maximum likelihood on the training fold's event times only, ignoring all covariates, so that every test subject receives an identical predicted survival probability this serves as a deliberate negative control, since a model with no subject-level discriminative information should achieve a concordance index of exactly 0.5; (iii) the BC-VAESM as described in Section 6, trained for 150 epochs per fold; and (iv) an ablation model, architecturally identical to the BC-VAESM in every respect (same encoder, same latent dimension d = 4, same training procedure and hyperparameters) but with its decoder modified to output the two parameters of a Weibull distribution in place of the three WEL parameters, isolating the specific contribution of the WEL parametric family. The evaluation metric is Harrell's concordance index (C-index), computed on the held-out fold using the predicted survival probability at the within-fold median event time as the risk score.

Implementation and reproducibility details All three neural models (Weibull-VAE, BC-VAESM, and the WEL-MLE negative control's underlying WEL fit) share the same encoder architecture: two fully connected layers (8 → 32 → 16 units, ReLU activations) mapping the eight standardized covariates (age, tumor size, tumor grade, positive lymph node count, progesterone receptor level, estrogen receptor level, hormonal therapy indicator, menopausal status) to a 4-dimensional latent Gaussian (μ, log σ²), and an identically structured decoder (latent + covariates → 32 → output) differing only in its output layer (3 units with softplus/sigmoid activations for (θ, α, c) in the BC-VAESM; 2 units with softplus activations for (shape, scale) in the Weibull-VAE ablation). Both are trained with the Adam optimiser (Kingma and Welling, 2014) at a fixed learning rate of 10⁻³, full-batch (the training fold, up to n = 549 observations, fits comfortably in memory), for a fixed 150 epochs with no early stopping (chosen based on the training-loss plateau already illustrated in Figure 7, Section 6, rather than a held-out validation criterion, since the five-fold outer loop already serves as the held-out evaluation), and gradient norm clipping at 5.0. Each of the five folds uses a distinct random seed (100 + fold for the Weibull-VAE, 200 + fold for the BC-VAESM, fold = 0, …, 4) so that results are exactly reproducible from the seed alone; the outer five-fold split itself uses scikit-learn's KFold with random_state=123. All experiments were run on CPU (no GPU was required, given the small covariate dimension and dataset size) using PyTorch 2.x. The Cox model uses lifelines's CoxPHFitter with a fixed ridge penalty of 0.1, chosen for numerical stability with correlated covariates rather than being tuned by cross-validation, since our purpose is a fair, standard baseline rather than an optimized competitor; we did not perform a hyperparameter search for the neural models, either using the single configuration above (chosen once, prior to seeing held-out results, by reference to the training dynamics of Section 6) for all five folds and all three neural/Cox models alike, precisely so that no model receives a comparative advantage from asymmetric tuning effort.

Table 13 and Figure 12 report the resulting five-fold cross-validated C-index, both as mean±standard deviation across folds and, in Table 14, as the individual per-fold values needed for the paired significance test below.

Table 13. Five-fold cross-validated C-index on GBSG2 (n = 686), from genuine model fits on real held-out data.

Table 14. Per-fold C-index values underlying Table 13, used for the significance tests in Table 15.

BC-VAESM attains the highest mean, though not significantly (Table 15).

Statistical significance Because all models are evaluated on the same five folds, the per-fold C-index values in Table 14 are a paired observation, and we test the BC-VAESM's improvement over each covariate-driven competitor with a paired t-test on these five paired differences (and, given the very small number of folds, also report the non-parametric paired Wilcoxon signed-rank test as a robustness check). Table 15 reports both. Neither comparison reaches conventional significance at n = 5 folds, and we state this plainly rather than letting the point estimates in Table 13 imply more than the data support: the evidence is directionally consistent with the WEL parametric family adding value within this architecture, but is not, by itself, statistically decisive at this sample size.

Table 15. Paired significance tests for the BC-VAESM's C-index improvement, from the fold-level differences in Table 14.

Interpretation Three findings emerge. First, the WEL-MLE negative control attains a C-index of exactly 0.500 in every fold, confirming both the correctness of our evaluation pipeline and the basic intuition that a model lacking subject-level covariate information cannot discriminate risk. Second, the BC-VAESM (mean C-index 0.684) modestly exceeds the penalized Cox model (mean 0.675); as Table 15 shows explicitly, this difference is not statistically significant at n = 5 folds (p = 0.29–0.31 depending on test), and we report it as a directionally promising but statistically inconclusive result at this sample size, rather than overstating its significance as the earlier version of this manuscript's mean-only comparison risked doing. Third, and most informative for assessing the genuine contribution of the WEL parametric family, the Weibull-decoder ablation (mean C-index 0.665) underperforms the full WEL-decoder BC-VAESM by 0.019 on average, across every fold but one, despite sharing an identical encoder architecture, latent dimension, and training procedure; this difference is also not statistically significant at n = 5 folds (p = 0.16–0.19). This comparison isolates the predictive value specifically attributable to the additional flexibility of the three-parameter WEL family relative to the two-parameter Weibull family within an otherwise unchanged deep generative architecture; we consider it suggestive, directionally consistent, empirical evidence for the practical utility of the BC-VAESM's parametric choice, and explicitly not a statistically confirmed one at the current sample size a larger, ideally multi-dataset or multi-cohort, evaluation would be needed to establish statistical significance, and we list this as a concrete direction for future work (Section 8).

Computational complexity and scalability Each BC-VAESM training run (150 epochs, full-batch, n ≈ 549 per fold) completed in under 10 seconds on a single CPU core in our environment; inference for a new subject is a single forward pass through the encoder and decoder, O(1) per subject and independent of training-set size, once trained, in contrast to the Cox model's partial-likelihood risk score, which in its basic form requires the training-set risk set at prediction time. Training cost scales linearly with the number of training subjects per epoch (a single forward/backward pass over the full batch) and is dominated by the incomplete gamma series evaluation (Section 6) inside the WEL log-likelihood, which we truncate at a fixed 60 terms per evaluation; for substantially larger cohorts than GBSG2's n = 686, mini-batch training (rather than the full-batch approach used here) would be a straightforward and standard modification, and the fixed-iteration series truncation means per-observation cost does not grow with dataset size. We have not tested the BC-VAESM beyond the dataset sizes examined in this article (n ≤ 686) and flag scalability to substantially larger cohorts (e.g. n > 10⁵) as untested rather than assumed.

Code and reproducibility The complete implementation (WEL distribution fitting, the differentiable incomplete-gamma function, the BC-VAESM and Weibull-VAE architectures and training loops, the Monte Carlo simulation harness, and the five-fold cross-validation pipeline) is implemented in Python (NumPy, SciPy, PyTorch, lifelines), and we will make the code publicly available in a repository at the time of publication, together with the exact random seeds used throughout this article, so that every table and figure in this manuscript can be independently regenerated.

8 Discussion

This article makes two contributions to lifetime data analysis, and we discuss both their strengths and their honestly-assessed limitations.

The first is the WEL distribution, a three-parameter lifetime model constructed via length-biased weighting of an Exponential-Lindley mixture. We initially conjectured, by analogy with composite-hazard constructions such as that of Ghitany et al. (2011), that this construction might yield multiple qualitatively distinct hazard shapes. Careful derivation of the mixture representation (Proposition 2.1) revealed instead that the resulting gamma-mixture structure shares a common scale parameter across both components, which by a classical result (Block and Savits, 1976) forces the hazard rate to be strictly increasing throughout the parameter space (Theorem 2.2). We view the identification and correction of this initially overstated structural property, in favor of the more modest but rigorously established IFR characterization, as itself a meaningful methodological contribution: the weight parameter c does not control hazard shape, but it does provide a genuine, independently verified, and empirically useful control over dispersion and tail weight, evidenced by the WEL's consistent log-likelihood improvement over the unweighted Lindley distribution, and its improvement over the two-parameter weighted Lindley distribution on GBSG2 and NCCTG Lung, across the four real datasets examined (Section 7.3). We do not extend this claim further than the evidence supports: on GBSG2 and on the guinea-pig data, the WEL's maximum likelihood estimate of α reaches the boundary of its parameter space, at which point the model's practical advantage over the plain weighted Lindley distribution effectively vanishes, and it is outperformed by the Generalized Gamma, Log-Normal, and/or Log-Logistic distributions.

The second contribution, the BC-VAESM, connects the WEL distribution to deep generative modeling by training a variational autoencoder whose decoder output is constrained to be a valid WEL parameter vector for every subject. The cross-validation results on GBSG2 show that this architecture is competitive with, and on the metric and dataset examined modestly exceeds a penalized Cox proportional hazards model, and that the specific choice of the WEL family (rather than a simpler Weibull family) contributes a measurable, though not large, predictive advantage within an otherwise identical architecture. We consider this ablation comparison, isolating the WEL family's contribution from the broader contribution of the variational autoencoder architecture itself, to be methodologically important and have prioritized including it over presenting a longer but less rigorous list of competing external baselines.

Limitations We identify the following limitations explicitly and without minimization. First, the WEL distribution cannot represent bathtub-shaped hazards, as demonstrated concretely by its poor fit to the AARSET dataset (Section 7.3); practitioners working with data known or suspected to exhibit non-monotonic hazard shapes should not use the WEL distribution. Second, even on datasets deliberately selected because their empirical hazard is increasing and therefore theoretically compatible with the WEL distribution's IFR property (the guinea-pig data, Section 7.1), the WEL distribution does not improve on the two-parameter Weibull, Gamma, or Log-Logistic distributions, or the three-parameter Generalized Gamma distribution, by AIC or BIC on any of the four datasets examined; its comparative advantage is specifically relative to other members of the Lindley mixture family (the unweighted and two-parameter weighted Lindley distributions), not relative to lifetime distributions in general. Third, on two of the four datasets (GBSG2 and the guinea-pig data) the maximum likelihood estimate of the mixing weight α reaches the boundary of its admissible range, indicating that α is only weakly identified once θ and c have been fitted, a practical estimability concern for the three-parameter model that applied users should be aware of. Fourth, the method of moments estimator exhibited a substantially reduced and configuration-dependent convergence rate (as low as 3–12% under Configuration III, and never exceeding 61% even at n = 300 under the configurations examined) in our simulation study, a practical limitation for this estimation approach that should be considered before recommending it for applied use. Fifth, the cross-validation evaluation of the BC-VAESM was conducted on a single dataset (GBSG2) with a sample size (n = 686) and event rate adequate for, but not large relative to, the number of folds and models compared; the modest improvement over Cox regression (mean C-index difference + 0.009, paired t-test p = 0.29, Wilcoxon p = 0.31; Table 15) is not statistically significant at this sample size and should be replicated, ideally with more folds or additional cohorts, before strong claims of superiority are made. Sixth, the BC-VAESM as implemented produces point predictions at inference time using the posterior mean of the latent representation rather than the full posterior averaging, for computational efficiency; full investigation of the calibration properties of the posterior predictive credible intervals described in Section 6.4 on additional datasets remains for future work.

Future work Natural extensions include a mixture of WEL components to accommodate multimodal or non-monotonic hazard behavior, thereby addressing the principal limitation identified above; a formal identifiability analysis of the mixing weight α to characterize the conditions under which its maximum likelihood estimate is well-separated from the boundary, motivated by the boundary estimates observed on GBSG2 and the guinea-pig data (Section 7.3); extension of the BC-VAESM to graph-structured or clustered survival data (e.g., patients sharing a treatment center); formal calibration assessment of the posterior predictive credible intervals via expected calibration error on larger cohorts; and application of the progressive Type-II censoring Bayesian machinery developed in Section 4.4 to data genuinely arising from a planned progressive censoring design (e.g., a reliability test with scheduled unit withdrawal), which remains distinct from the random right-censoring already incorporated into the GBSG2 and NCCTG Lung analyses of Section 7 and was beyond the scope of the present study.

9 Conclusion

We have introduced the Weighted Exponential-Lindley (WEL) distribution, a three-parameter continuous lifetime model obtained by applying a length-biased weight function to a mixture of exponential and Lindley components sharing a common scale parameter. We proved that this construction yields a distribution with a strictly increasing hazard rate throughout its parameter space, converging to the scale parameter as x → ∞, and demonstrated that the additional weight parameter provides an independent, continuously tunable control over dispersion, skewness, and kurtosis that is unavailable in the unweighted Lindley distribution or the existing two-parameter weighted Lindley distribution. Comprehensive distributional properties were derived in closed or computationally tractable form, and maximum likelihood, method-of-moments, and Bayesian estimation under progressive Type-II censoring was developed and compared through a genuine Monte Carlo simulation study, which showed a substantial and consistent advantage for Bayesian estimation of the mixing-weight parameter specifically.

The Bayesian-Calibrated Variational Autoencoder Survival Model (BC-VAESM) connects this distributional theory to deep generative modeling by training an encoder-decoder architecture whose decoder output is constrained to be a valid WEL parameter vector for every subject. Five-fold cross-validation on the real GBSG2 breast cancer dataset showed that the BC-VAESM is competitive with a penalized Cox proportional hazards model, and an architecturally controlled ablation isolated a genuine, measurable contribution of the WEL parametric family relative to a simpler Weibull alternative within the same deep learning architecture.

Balancing strengths against limitations, as we believe, a conclusion should: on the strength side, the WEL distribution is analytically tractable (closed-form density, survival function, and moments via the gamma-mixture representation), its hazard-shape property is proved rigorously rather than assumed, its parameters have direct practical interpretations (Section 7.3), and it measurably improves on the unweighted Lindley distribution, and on the weighted Lindley distribution on two of four real datasets, across every dataset examined; the BC-VAESM adds calibrated latent-variable uncertainty quantification to an otherwise standard covariate-driven survival architecture at a modest, characterized computational cost (Section 7). On the limitation side, we state these without minimization: the WEL distribution's failure to improve on the Weibull, Gamma, Log-Normal, Log-Logistic, Birnbaum-Saunders, or Generalized Gamma distributions by AIC or BIC on any of the four real datasets examined, its inability to fit the bathtub-shaped AARSET data, the boundary (weakly identified), the maximum likelihood estimate of the mixing weight α observed on GBSG2 and on a dataset (guinea-pig survival times) selected specifically to match its theoretical hazard shape, and the modest, not statistically significant (p > 0.16 by either a paired t-test or the Wilcoxon signed-rank test, Table 15) margin by which the BC-VAESM exceeds Cox regression and the Weibull-decoder ablation are all real and are all reported above exactly as obtained.

Practical implications and limitations for applied survival analysis For a practitioner deciding whether to use the WEL distribution or the BC-VAESM on a new dataset, our results suggest the following, stated as guidance rather than a guarantee: the WEL distribution is a reasonable default when (i) the data are known or suspected to have an increasing hazard (it should not be used, and will not fit well, if a bathtub or decreasing hazard is plausible, Section 7.3), and (ii) a simple, closed-form, readily interpretable parametric model is preferred over a possibly better-fitting but less interpretable three-parameter alternative such as the Generalized Gamma; where maximal goodness-of-fit is the sole objective, our own results show that trying several two- and three-parameter competitors and comparing by AIC/BIC, exactly as we do throughout Section 7.3, remains good practice, since no single model (WEL included) dominated across all four datasets, we examined. The BC-VAESM is most useful when an individual (not merely population-average) predictive uncertainty is required, and a moderately sized, richly covariate-annotated dataset (order hundreds to a low thousands of subjects, based on the one dataset examined here) is available for training; we have not evaluated it on datasets substantially larger or smaller than GBSG2's n = 686, on datasets with markedly different censoring proportions, or on data exhibiting a non-monotonic hazard, and recommend against extrapolating our results to those settings without further validation.

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: Ayman Baklizi, Yarmouk University Faculty of Sciences, Jordan

Reviewed by: Ibrahim E. Ragab, Alexandria University, Egypt

Reza Azimi, University of Mohaghegh Ardabili, Iran