Bounded Data Modeling with the Extended Bradford Distribution: Modal Regression Approach and Applications

1  Introduction

Regression analysis is a fundamental tool for examining the relationship between a set of explanatory variables and a response variable. The standard regression approaches focus on modeling the conditional expectation of the response variable, representing the central tendency of its association with the covariates. However, relying on the mean structure may be insufficient in practice, as such methods can lack robustness, particularly in datasets affected by outliers or characterized by heavy-tailed distributions (see [1,2]).

Models based on the conditional mode can provide a more useful representation of the underlying relationship than methods based on the conditional mean or median when the conditional distribution of the response variable given the covariates displays skewness or heavy tails. Reference [3] introduced two parametric modal regression models for bounded responses based on the beta distribution and provided diagnostic tools. Reference [4] presented a flexible parametric modal regression framework based on a heterogeneous Gumbel mixture distribution. Reference [5] developed a Bayesian modal regression using an unimodal distribution capable of accommodating symmetry, asymmetry, and varying tail behaviours through its shape and scale parameters. Reference [6] applied modal regression to neuroimaging data, demonstrating that traditional mean regression failed to capture meaningful associations due to the skewed and heavy-tailed nature of the cognitive assessment data. Reference [7] proposed three new modal regression models for bounded data based on unit-Gamma, kumaraswamy, and unit-Gompertz distributions, and showed that these models offer flexible and effective alternatives for modeling bounded responses under different data characteristics.

The primary motivation of this study arises from the increasing need for robust regression frameworks capable of effectively modeling bounded response variables in the presence of skewness, heavy tails, and outliers. While modal regression has recently gained attention as a robust alternative to mean-based approaches, its practical implementation remains limited due to the lack of suitable parametric distributions with analytically tractable modes.

Although several bounded distributions exist in the literature, most of them do not provide closed-form expressions for the mode, which restricts their applicability in a modal regression model (see [7]). Some of the recently proposed bounded distributions are examined. For instance, although the distribution of [8] is flexible, its mode cannot be obtained in closed form. The Poisson-unit-Weibull distribution of [9] is suitable for quantile regression, but it cannot be used for conditional mean and mode modeling of the response variable. Although the unit inverse exponentiated Lomax distribution of [10] is quite satisfactory at modeling various hazard forms, it lacks closed-form equations for mean and mode values. So, in the literature, many flexible distributions are available; however, their use in generalized linear models, time series, or engineering modeling is limited. Recently, Reference [11] introduced a one-parameter bounded distribution, called the Bradford distribution, with the following probability density function (pdf):

f(y)=θlog⁡(θ+1)(θy+1),(1)

where 0<y<1 and θ>0 is the scale parameter. While the Bradford distribution is notable for its simple, univariate structure, its use in modeling real-world datasets is quite limited. The Bradford distribution is right-skewed and has a monotonically increasing hazard rate function (hrf). In contrast, real-world datasets frequently exhibit different types of hrfs and left-skewness. To increase the flexibility of the Bradford distribution and improve its ability to model datasets with different characteristic structures, the extended Bradford (EB) distribution is proposed. The EB distribution can be used to model data that is skewed in both directions. It is also capable of modeling a variety of hrfs such increasing and bathtub. The EB distribution is suitable for modal regression models since its mode value has an explicit mathematical expression.

The proposed distribution differs from existing generalized bounded distributions in several important aspects. First, it preserves the analytical simplicity of the classical Bradford distribution while introducing an additional shape parameter that substantially increases flexibility. Second, the model is capable of modeling both left and right-skewed pdfs together with various hazard rate structures. Third, the mode of the distribution admits an explicit analytical expression, which enables a natural and computationally convenient reparameterization for modal regression purposes. These properties provide important practical advantages over several existing bounded models, including the beta and Kumaraswamy-based modal regression models. Therefore, the proposed approach offers flexibility, interpretability, and computational simplicity.

The presented study makes several key contributions to the literature. It introduces a new two-parameter bounded distribution that generalizes the classical Bradford distribution. It points out the key mathematical properties of the proposed distribution, including its moments, hazard rate behavior, and structural relationship with the beta distribution. It provides a closed-form expression for the mode, making the distribution particularly suitable for modal regression modeling. It develops a modal regression framework based on the EB distribution, including parameter estimation via maximum likelihood and diagnostic tools based on randomized quantile residuals. It implements the proposed model within an interactive R Shiny application, enhancing its accessibility and practical usability.

The rest of this paper is structured in the following way: In Section 2, the EB distribution is introduced, and its fundamental properties are presented. Section 3 is devoted to parameter estimation, where the maximum likelihood estimation procedure and its asymptotic properties are discussed along with a simulation study. In Section 4, the modal regression framework based on the EB distribution is developed, and the performance of the estimators is evaluated through Monte Carlo simulations. Section 5 illustrates the practical applicability of the model using real data examples, including both univariate and regression settings. In Section 6, an interactive implementation of the proposed method via an R Shiny application is described. Section 7 contains the limitations and possible future research.

2  Extended Bradford Distribution

Applying X=Y1α transformation where Y denotes a random variable with the pdf in (1), the pdf of the EB distribution is given by

f(x)=αθxα−1log⁡(θ+1)(θxα+1), 0<x<1,(2)

where α>0 is the shape parameter. The EB distribution reduces to the Bradford distribution for α=1. The cumulative distribution function (cdf) of the EB distribution is

F(x)=log⁡(1+θxα)log⁡(1+θ).(3)

The quantile function (qf) of the EB distribution is

Q(p)=((1+θ)p−1θ)1α,(4)

where 0<p<1. The qf of the EB distribution is obtained in closed form. Therefore, the inverse transform method is easily implemented to generate random variables from the EB distribution. Fig. 1 shows the pdf shapes of the EB distribution according to different parameter values. The EB distribution stands out as a distribution that can be used to model extremely left or right-skewed data. The hrf of X is

h(x)=αθxα−1(log⁡(1+θ)−log⁡(1+θxα))(1+θxα).(5)

Figure 1: Pdf shapes of the EB distribution.

The hrf shapes of the EB distribution are displayed in Fig. 2. The distribution has two failure shapes: bathtub and increasing.

Figure 2: Hrf shapes of the EB distribution.

In the following subsection, the properties of the EB distribution are discussed in detail.

2.1 Connection with the Beta Distribution

The EB distribution has a structural connection with the beta distribution. The pdf of the beta distribution for β=1 is

g(x)=αxα−1,0<x<1.(6)

The pdf of the proposed model can be rewritten as

f(x)=g(x)θlog⁡(1+θ)(θxα+1),(7)

which reveals that the EB distribution can be interpreted as a weighted version of the Beta(α,1) distribution. In this representation, the weighting function is

w(x)=1θxα+1,(8)

and log⁡(1+θ)/θ is the normalizing constant. This representation shows that the proposed distribution belongs to the class of weighted distributions, where the baseline density is the Beta(α,1) distribution and the weight function introduces additional flexibility in modeling skewness and tail behavior. Furthermore, as θ→0,

f(x)→αxα−1,(9)

which is the pdf of the Beta(α,1) distribution.

2.2 Properties

Proposition 1. The rth raw moment of X is

E(Xr)=θαlog⁡(1+θ)(α+r)F12(1,1+rα;2+rα;−θ),(10)

where F12(⋅) denotes the Gauss hypergeometric function.

Proof. The r-th raw moment is

E(Xr)=∫01xrf(x;α,θ)dx.(11)

Substituting the pdf in (11), we get

E(Xr)=αθlog⁡(1+θ)∫01xr+α−11+θxαdx.(12)

Let u=xα, du=αxα−1dx and xα−1dx=1αdu. Also, xr+α−1dx=urα1αdu. Since x∈(0,1), u∈(0,1). Applying these substitution in (12), we have

E(Xr)=θlog⁡(1+θ)∫01urα1+θudu.(13)

Let rα=a−1, the integral in (13) is

∫01ua−11+θu.du.(14)

Using the integral representation of the Gauss hypergeometric function, we obtain

∫01ua−11+θudu=1aF12(1,a;a+1;−θ).(15)

Applying (15) in (13), the Eq. (13) becomes

E(Xr)=θlog⁡(1+θ)α F12(1,a;a+1;−θ).(16)

Substituting a=1+rα in (16), we get

E(Xr)=θαlog⁡(1+θ)(α+r)F12(1,1+rα;2+rα;−θ).(17)

□

Proposition 2. The mean of X is

E(X)=θ−1/αlog⁡(1+θ)Bθ1+θ(1+1α,−1α),(18)

where Bz(a,b) denotes the incomplete beta function, given by [12].

Proof. The mean of X is

E(X)=∫01xf(x)dx=αθlog⁡(1+θ)∫01xα1+θxαdx.(19)

Two transformations are used for the integration part of (19). First, let u=xα,du=αxα−1dx. Then, the integral becomes

E(X)=θlog⁡(1+θ)∫01u1/α1+θudu.(20)

For the second transformation, let y=θu1+θu, u=yθ(1−y) and du=1θ(1−y)2dy. When u=0, y=0; when u=1, y=θ1+θ. Substituting these change of variables in (20), we have

E(X)=θ−1/αlog⁡(1+θ)∫0θ1+θy1/α(1−y)−1/α−1dy.(21)

The integral representation of the incomplete beta function is

Bx(a,b)=∫0xta−1(1−t)b−1dt.(22)

It is clear that the integration part in (21) is the incomplete beta function with the following parameters x=θ/(θ+1), a=1/α+1 and b=−1/α.□

The expressions obtained in Propositions 1 and 2 are closely related. In Proposition 1, the rth raw moment of the EB distribution is expressed in terms of the Gauss hypergeometric function. In particular, setting r=1 yields the mean of the distribution. On the other hand, in Proposition 2, the mean is expressed in terms of the incomplete beta function. Applying y=θu/(1+θu) transformation in (15), we obtain

∫01ua−11+θudu=θ−a∫0θ1+θya−1(1−y)−ady=θ−aBθ1+θ(a,1−a),(23)

which is the incomplete beta function representation of (15). Therefore, the expressions given in Propositions 1 and 2 are mathematically consistent. Using the similar approach in Proposition 2, the second raw moment of X is

E(X2)=θ−2/αlog⁡(1+θ)Bθ1+θ(1+2α,−2α).(24)

So, the variance of X is

Var⁡(X)=θ−2/αlog⁡(1+θ)Bθ1+θ(1+2α,−2α)−[θ−1/αlog⁡(1+θ)Bθ1+θ(1+1α,−1α)]2.(25)

Proposition 3. The moment generating function (mgf) of X is

MX(t)=θlog⁡(1+θ)∑k=0∞tkk!11+kαF12(1,1+kα;2+kα;−θ).(26)

Proof. The mgf is defined as

MX(t)=αθlog⁡(1+θ)∫01xα−1etx1+θxαdx.(27)

Using the power series expansion for exp⁡(tx)=∑k=0∞tkxkk!, the mgf is rewritten as

MX(t)=αθlog⁡(1+θ)∑k=0∞tkk!∫01xk+α−11+θxαdx.(28)

Using u=xα transformation, we obtain

MX(t)=θlog⁡(1+θ)∑k=0∞tkk!∫01ukα1+θudu.(29)

As in Proposition 1, the integral representation of the Gauss hypergeometric function can be used. Let a=1+kα, the integration part in (29) can be written as

∫01ua−11+θudu=1aF12(1,a;a+1;−θ).(30)

Inserting (30) in (29), we get

MX(t)=θlog⁡(1+θ)∑k=0∞tkk!1aF12(1,a;a+1;−θ).(31)

Replacing a=1+kα in (31), the mgf is

MX(t)=θlog⁡(1+θ)∑k=0∞tkk!11+kαF12(1,1+kα;2+kα;−θ).(32)

□

Proposition 4. If α≤1, the pdf is strictly decreasing on (0,1) and the mode is at x=0. If α>1, the pdf is unimodal and the mode is

x∗=(α−1θ)1α,(33)

where (α−1)/θ<1. Otherwise the mode occurs at x=1.

Proof. Since α>0 and θ>0, the normalizing constant αθlog⁡(1+θ) does not affect the mode. Hence it is sufficient to maximize the below function to determine the mode of the EB distribution.

g(x)=xα−11+θxα.(34)

The derivative of (34) is

g′(x)=xα−2[(α−1)−θxα](1+θxα)2.(35)

Since the denominator is strictly positive and xα−2>0, the mode of the EB distribution is obtained by solving (α−1)−θxα=0. So, the mode is

x∗=(α−1θ)1α.(36)

If α≤1, then (α−1)≤0, so the equation has no positive solution in (0,1). In this case, g′(x)<0 for all x∈(0,1) and the pdf is strictly decreasing; hence the mode is at x=0.

If α>1, then (α−1)/θ>0, so a unique positive critical point exists. If (α+1)/θ<1, then x∗∈(0,1) and this point yields the unique maximum since g′(x) changes sign from positive to negative. If α−1θ≥1, then the critical point lies outside the support, and g′(x)>0 throughout (0,1), so the pdf is increasing and the mode occurs at x=1.

□

Proposition 5. The EB distribution is log-convex for α≤1.

Proof. Consider the log-pdf of the EB distribution

log⁡f(x)=log⁡α+log⁡θ−log⁡log⁡(1+θ)+(α−1)log⁡x−log⁡(1+θxα).(37)

The second derivative of (37) with respect to x is

d2dx2log⁡f(x)=−α−1x2−θαxα−2[(α−1)−θxα](1+θxα)2.(38)

For α≤1, we have

−α−1x2≥0,(α−1)−θxα<0∀x∈(0,1).(39)

Hence the second term is positive, and we have

d2dx2log⁡f(x)>0,∀x∈(0,1).(40)

Thus the log-pdf is convex, and the pdf is log-convex.

For α>1, the expression (α−1)−θxα is positive for sufficiently small x and becomes negative for sufficiently large x. Therefore, the second derivative changes sign on (0,1). Consequently, the log-pdf is neither concave nor convex on (0,1).

□

Proposition 6. The following results are obtained for the shape of the hrf:

(i)   If 0<α<1, then h(x) has a bathtub shape; that is, it is initially decreasing, attains a minimum, and then increases to infinity.

(ii)   If α=1, then h(x) is increasing on (0,1).

(iii)   If α>1, then h(x) has a J-shape that starts from zero and goes to infinity.

Proof. The hrf is given by

h(x)=αθxα−1(1+θxα)[log⁡(1+θ)−log⁡(1+θxα)].(41)

When x→0+ and x→1−, the limits of the hrf are given by

limx→0+h(x)={∞,0<α<1,θlog⁡(1+θ),α=1,0,α>1,(42)

and

limx→1−h(x)=∞.(43)

Consider the derivative of log⁡h(x), given by

ddxlog⁡h(x)=α−1x−θαxα−11+θxα[1+1log⁡(1+θ)−log⁡(1+θxα)].(44)

The last term of (44) is strictly greater than 1 for all x∈(0,1). The three cases are examined. The first case is for 0<α<1. As x→0+, we have

α−1x→−∞.

So, we derive

ddxlog⁡h(x)<0,

and h(x) is decreasing near zero. As x→1−, we have

ddxlog⁡h(x)→−∞,

and h(x)→∞. Since h(x) first decreases and goes infinity, it has a bathtub shape 0<α<1.

The second case is for α=1. The hrf is

h(x)=θ(1+θx)[log⁡(1+θ)−log⁡(1+θx)].(45)

Taking the first derivative of (45), we get

ddxh(x)>0for all x∈(0,1).

Hence, h(x) is strictly increasing for α=1. The third case is for α>1. As x→0+, we have h(x)→0 and

ddxlog⁡h(x)∼α−1x>0.

So, h(x) increases near zero. As x→1−, h(x)→∞. Therefore, h(x) increases from zero to infinity. So, h(x) has J-shaped for α>1.

□

Proposition 7. The Rényi entropy of X is

Rδ(α,θ)=11−δlog⁡[1α(αθlog⁡(1+θ))δ1aF12(δ,a;a+1;−θ)],(46)

where

a=δ(α−1)+1α.

Proof. The Rényi entropy is

Rδ(α,θ)=11−δlog⁡(∫01f(x)δdx).(47)

For the EB distribution, the integral part is

∫01f(x)δdx=(αθlog⁡(1+θ))δ∫01xδ(α−1)(1+θxα)−δdx.(48)

Using u=xα transformation, the integral is rewritten as follows:

∫01f(x)δdx=1α(αθlog⁡(1+θ))δ∫01uδ(α−1)+1α−1(1+θu)−δdu.(49)

As in Proposition 1, the integral in (49) can be expressed in terms of the Gauss hypergeometric function, as follows:

∫01ua−1(1+θu)−δdu=1aF12(δ,a;a+1;−θ),(50)

where

a=δ(α−1)+1α.(51)

So, the Rényi entropy of X is

Rδ(α,θ)=11−δlog⁡[1α(αθlog⁡(1θ))δ1aF12(δ,a;a+1;−θ)].(52)

□

For a detailed shape analysis of the EB distribution, its skewness and kurtosis values are obtained via computational implementation of the Gauss hypergeometric function in the gsl package of R. The computational results in Fig. 3 indicate that an increase in α leads to a decrease in both skewness and kurtosis, whereas an increase in θ yields the opposite effect.

Figure 3: Skewness and kurtosis values of the EB distribution.

3  Estimation

The parameter estimation for the EB distribution via the maximum likelihood estimation (MLE) method is detailed in this section. Furthermore, the asymptotic efficiency of the ML estimators is evaluated through a comprehensive simulation study.

3.1 Maximum Likelihood

The log-likelihood function of the EB distribution is

ℓ(α,θ)=nlog⁡α+nlog⁡θ+(α−1)∑i=1nlog⁡xi−nlog⁡log⁡(1+θ)−∑i=1nlog⁡(1+θxiα).(53)

Differentiating (53) with respect to α and θ, we obtain

∂ℓ∂α=nα+∑i=1nlog⁡xi−∑i=1nθxiαlog⁡xi1+θxiα,(54)

∂ℓ∂θ=nθ−n(1+θ)log⁡(1+θ)−∑i=1nxiα1+θxiα.(55)

The ML estimators (α^,θ^) are obtained by solving the following nonlinear system:

∂ℓ∂α=0,∂ℓ∂θ=0.

Since these equations involve nonlinear terms, the closed-form solutions do not exist. Therefore, iterative optimization algorithms should be employed. The elements of the observed information matrix are given by

I(α,α)=∂2ℓ∂α2=−nα2−∑i=1nθxiα(log⁡xi)2(1+θxiα)2,(56)

I(θ,θ)=∂2ℓ∂θ2=−nθ2+n[(1+θ)log⁡(1+θ)−1](1+θ)2log2⁡(1+θ)−∑i=1nxi2α(1+θxiα)2,(57)

I(α,θ)=∂2ℓ∂α∂θ=−∑i=1nxiαlog⁡xi(1+θxiα)2.(58)

When the standard regularity conditions hold, the MLEs for the EB distribution are consistent and asymptotically normal. Let Θ={(α,θ):α>0,θ>0} be the parameter space. Given that the support of the distribution is parameter-independent and its log-likelihood function is twice continuously differentiable with respect to (α,θ), the regularity conditions are satisfied. Moreover, model identifiability is ensured since distinct parameter vectors result in distinct PDFs. Therefore, based on classical likelihood theory, it follows that

n(α^−αθ^−θ)→dN2(0,I(α,θ)−1).

3.2 Simulation

The performance of the MLEs in terms of bias, mean squared error (MSE), and mean relative error (MRE) is discussed under different sample sizes and parameter vectors. The simulation is repeated 1000 times. The used parameter vectors are (α=2,θ=2), (α=0.5,θ=0.5), and (α=1.5,θ=0.5).

The simulation results given in Table 1 reveal clear and consistent results for all parameter settings. As the sample size increases from n=100 to n=5000, all performance measures improve substantially. Specifically, the bias of both α^ and θ^ decreases monotonically toward zero, indicating that the estimators are asymptotically unbiased. Similarly, the MSE values decrease significantly with increasing sample size, confirming the consistency of the MLEs. Moreover, the MRE values approach 1 as n increases, suggesting that the estimators become increasingly accurate in relative terms.

A notable difference is observed between the estimation performance of α and θ. The estimator of α performs considerably better across all scenarios, exhibiting smaller bias, lower MSE, and MRE values very close to 1 even for relatively small sample sizes. In contrast, the estimator of θ shows substantially higher bias and MSE, particularly when n=100. Additionally, the MRE values for θ are significantly larger than 1 in small samples, indicating poor relative accuracy. This suggests that θ is more difficult to estimate and may be more sensitive to sampling variability. Therefore, bias-reduction techniques or alternative estimation methods should be considered.

Note that all computations are conducted in the R software. The MLEs are obtained using the optim function with the Nelder-Mead optimization algorithm. To improve numerical stability and reduce sensitivity to initial values, multiple starting values are considered in the optimization procedure. The observed information matrix is numerically evaluated at the MLEs of the EB distribution and used to compute standard errors and asymptotic confidence intervals. Throughout the simulation experiments, the convergence is achieved in almost all replications.

4  Modal Regression

The mode of the EB distribution is (for α>1)

x∗=(α−1θ)1α.(59)

The EB distribution is re-parametrized by means of the its mode. Let θ=(α−1)μ−α, the pdf of the mode-parametrized EB distribution is

f(y)=α(α−1)μ−αyα−1log⁡((α−1)μ−α+1)((α−1)μ−αyα+1).(60)

The pdf in (60) is denoted as EB(μ,α), where α>0, μ∈(0,1), and mode(Y)=μ. The shapes of the pdf in (60) are displayed in Fig. 4.

Figure 4: Pdf shapes of the mode-parametrized EB distribution.

Let yi,i=1,2,...,n be the realizations of the random variable with the pdf in (60). The conditional mode of the parametrized EB distribution is modeled with the logit link function

μi=exp⁡(xiγT)1+exp⁡(xiγT),(61)

where γ=(γ0,γ1,γ2,…,γp)T is unknown regression parameter vector and xi=(1,xi1,xi2,xi3,…,xip)T is the ith vector of known covariates. Since the mode of the EB distribution is within the range of (0, 1), the logit link function is used.

Substituting the link function (61) in (60), the log-likelihood function of the EB modal regression is

ℓ(γ,α)=nlog⁡(α(α−1))−α∑i=1nlog⁡(μi)+(α−1)∑i=1nlog⁡(yi)−∑i=1nlog⁡(log⁡((α−1)μi−α+1))+∑i=1nlog⁡((α−1)μi−αyiα+1),(62)

where μi is defined in (61). As expected, there is no explicit solution for the ML estimators of (γ,α). For this reason, direct maximization of (62) is required to obtain (γ^,α^). There are many available packages in several software to maximize the likelihood function. Here, the optim function in R is preferred. Based on the observed information matrix, the standard errors and asymptotic confidence intervals of the parameters are constructed. The model accuracy is evaluated by the randomized quantile residual (see [13])

ri=Φ−1(pi),

where pi=Fμ^i,α^(yi) and Φ−1 is the qf of the N(0,1) distribution. When the fitted model adequately captures the data, ri follows the N(0,1) distribution.

4.1 Simulation of EB Modal Regression Model

The asymptotic behavior of the ML estimates for the EB modal regression model is investigated through a Monte Carlo simulation study. The dependent variable is generated using the logit link function defined as

logit(μi)=γ0+γ1xi1+γ2xi2.(63)

Note that μi in (63) can be rewritten as in (61). Two different parameter settings are considered:

•   Case I: γ0=−2,γ1=−1,γ2=0.5,α=2,

•   Case II: γ0=−3,γ1=1.5,γ2=−1.5,α=2.

The covariates xi1 and xi2 are independently generated from the standard uniform distribution. For each configuration, the simulation is replicated 1000 times. The performance of the ML estimators is evaluated in terms of the empirical mean, bias, and MSE, and the results are reported in Table 2.

The results indicate that the ML estimators perform satisfactorily in both scenarios. As the sample size increases, the empirical means of the estimates approach the true parameter values, demonstrating consistency. The bias values are relatively small even for moderate sample sizes and tend to decrease as n increases.

Moreover, the MSE values exhibit a clear decreasing pattern with increasing sample size, confirming the efficiency of the estimators. This behavior is observed consistently across all parameters in both Case I and Case II. In particular, the estimation of the shape parameter α shows relatively larger bias and MSE for smaller sample sizes; however, its performance improves significantly as the sample size grows. In summary, these numerical findings offer robust empirical support for the consistency and asymptotic efficiency of the parameter estimates in the EB modal regression model.

5  Applications

5.1 Univariate Case

Reference [14] introduced the log-cosine-power (LCP) distribution and demonstrated its flexibility using a real data set on failure times. This data set, originally analyzed in [14], has also been considered in comparison with several competing unit distributions, including the unit Teissier distribution [15], transmuted unit Rayleigh distribution [16], Topp-Leone distribution [17], unit-exponential distribution [18], and unit Burr XII distribution [19]. According to [14], the LCP distribution provides the best fit among these models, as evidenced by the lowest Akaike information criterion (AIC) and Bayesian information criterion (BIC) values, along with the highest p-value of the Kolmogorov-Smirnov (KS) test. Specifically, the reported AIC, BIC, and p-value for the LCP model are −7.0748, −3.1721, and 0.5635, respectively.

In this study, the same data set is reanalyzed to assess the performance of the proposed EB distribution in (2), and compared with well-known alternative models. The competing models considered are the beta and Kumaraswamy (Kum) [20] models. The estimation results and goodness-of-fit statistics are presented in Table 3.

The results indicate that the EB model outperforms the competing models. In particular, it yields the smallest AIC and BIC values. Furthermore, the EB model produces the lowest KS statistic and the highest corresponding p-value, providing strong evidence that it offers the closest agreement with the empirical distribution of the data.

The graphical analysis presented in Fig. 5 provides additional support in favour of the EB distribution. The fitted PDF of the EB distribution is close to the empirical histogram and successfully captures both the central tendency and the tail behavior. Overall, both the numerical and graphical findings provide adequate evidence that the EB distribution offers a more flexible and accurate modeling for the considered failure time data.

Figure 5: Fitted pdfs of the considered models.

5.2 Regression with Covariates

The water quality (WQ) indicator of the OECD (Organisation for Economic Co-operation and Development) countries is modelled with the air pollution (AP) and long-term unemployment rate (LTUR) indicators. The EB modal regression model is compared with the beta and gamma modal regression models. The results are reported in Table 4.

The coefficient of AP is negative and statistically significant for all fitted models. It indicates that higher levels of AP are associated with lower levels of WQ. In other words, an increase in AP leads to a deterioration in WQ, which is consistent with environmental theory and empirical expectations. Similarly, the LTUR exhibits a negative and significant effect on WQ in all models. It suggests that higher unemployment rates are linked to poorer environmental outcomes, possibly reflecting reduced public and private investment in environmental protection. Overall, the estimated coefficients provide consistent evidence that both environmental and socio-economic factors play a crucial role in determining WQ for OECD countries.

Additionally, it can be seen that the EB modal regression model has lower AIC and BIC values than the other two models. Therefore, the EB model yields better results than the beta and gamma modal regression models for this dataset.

The results of the KS test, applied to the randomized quantile residuals of the fitted regression models, are presented in Table 5. The residuals obtained from all fitted models follow a standard normal distribution. However, the p-value obtained for the EB model is higher than that of the other two models.

Fig. 6 shows the quantile-quantile (QQ) plot for the residuals. It is worth noting that the theoretical quantile values in the QQ plot obtained for the EB model are closer to the empirical values than those for the other two models.

Figure 6: QQ plots of the residuals.

6  Implementation of EB Modal Regression Model via R Shiny

An interactive web-based application is developed using the R Shiny framework to facilitate the implementation of the proposed model. The application, available at https://gazistat.shinyapps.io/EBMR/, is designed to provide an accessible and user-friendly environment for performing estimation, model evaluation, and diagnostic analysis of the EB model.

The interface is structured into four main panels: Data, Regression, Diagnostics, and Univariate Fit. Each panel is designed to guide the user through a specific stage of the analysis.

The Data panel allows users to upload datasets in commonly used formats such as .csv and .xlsx (see Fig. 7). Within this panel, users can select the dependent and independent variables and identify categorical covariates. For categorical variables, the application enables the specification of reference levels, ensuring proper model interpretation. In addition, a validation mechanism checks whether the dependent variable lies within the unit interval and provides a warning if the condition is violated.

Figure 7: Data upload and variable selection panel.

The Regression panel presents the results obtained from the estimation procedure (see Fig. 8). It reports parameter estimates along with their standard errors, test statistics, and significance levels. Furthermore, commonly used model selection criteria is provided to support model comparison and evaluation.

Figure 8: Regression panel.

The Diagnostics panel is devoted to assessing model adequacy (see Fig. 9). It provides graphical and numerical tools based on randomized quantile residuals. In particular, the QQ plot is displayed to visually evaluate the agreement between empirical and theoretical distributions. In addition, the KS test is reported to assess the goodness-of-fit. The panel also allows users to export the diagnostic plot in PDF and EPS formats for reporting purposes.

Figure 9: Diagnostics panel.

The Univariate Fit panel enables users to fit the proposed distribution to univariate data (see Fig. 10). Users can upload a dataset, select a variable, and specify initial parameter values. The panel reports parameter estimates, information criteria, and goodness-of-fit measures. It also provides graphical outputs, including fitted density curves, hrfs, survival functions, and probability plots, which help to visually assess the adequacy of the fitted distribution.

Figure 10: Univariate data fitting panel.

7  Conclusion and Limitations

In this study, a new bounded distribution, referred to as the EB distribution, is introduced as a flexible extension of the classical Bradford distribution. The proposed distribution incorporates an additional shape parameter, allowing it to capture a wide range of distributional behaviors, including various forms of skewness and hazard rate structures. Its structural properties, as well as its connection with the beta distribution, make it both theoretically and practically convenient. A key advantage of the EB distribution lies in the availability of an explicit expression for the mode, which makes it suitable for modal regression modeling. Building on this property, a modal regression framework is developed, enabling the modeling of conditional modes through a parametric and interpretable structure. The estimation procedures based on maximum likelihood are shown to be consistent and efficient as the sample size increases. The practical utility of the EB distribution and its associated regression model is demonstrated through real data applications. The results indicate that the proposed model provides improved fit and greater flexibility compared to several well-known alternatives. In addition, the implementation of the model through an interactive R Shiny application enhances its accessibility in applied research.

Although the EB distribution is flexible and analytically tractable, the distribution has some limitations. First, the numerical optimization required for MLEs may become computationally difficult for large-scale applications. The proposed model is currently developed for observed bounded response variables and does not directly accommodate censored or truncated data structures that are common in survival and reliability studies. Additionally, the modal regression is studied only for low-dimensional settings, and extensions to high-dimensional covariate spaces may require regularization or penalized estimation techniques. Overall, the EB distribution offers a valuable contribution to the class of bounded distributions and modal regression methodology. Future research may focus on extending the proposed framework to more complex settings such as variable dispersion models, comparison of alternative link functions, or robust estimation techniques.

Acknowledgement: Not applicable.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Emrah Altun, Christophe Chesneau, Atacan Erdis; methodology, Emrah Altun, Christophe Chesneau, Atacan Erdis; software, Emrah Altun, Christophe Chesneau, Atacan Erdis; validation, Emrah Altun, Christophe Chesneau, Atacan Erdis; investigation, Emrah Altun, Christophe Chesneau, Atacan Erdis; writing—original draft preparation, Emrah Altun, Christophe Chesneau, Atacan Erdis; writing—review and editing, Emrah Altun, Christophe Chesneau, Atacan Erdis; visualization, Emrah Altun, Christophe Chesneau, Atacan Erdis. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data supporting the findings of this study are available from the corresponding author upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Sager TW, Thisted RA. Maximum likelihood estimation of isotonic modal regression. Ann Statist. 1982;10(3):690–707. doi:10.1214/aos/1176345865. [Google Scholar] [CrossRef]

2. Collomb G, Hardle W, Hassani S. A note on prediction via estimation of the conditional mode function. J Stat Plan Inference. 1986;15:227–36. doi:10.1016/0378-3758(86)90099-6. [Google Scholar] [CrossRef]

3. Zhou H, Huang X, Alzheimer’s Disease Neuroimaging Initiative. Parametric mode regression for bounded responses. Biom J. 2020;62(7):1791–809. doi:10.1002/bimj.202000039. [Google Scholar] [PubMed] [CrossRef]

4. Liu Q, Huang X, Zhou H. The flexible gumbel distribution: a new model for inference about the mode. Stats. 2024;7(1):317–32. [Google Scholar]

5. Liu Q, Huang X, Bai R. Bayesian modal regression based on mixture distributions. Comput Stat Data Anal. 2024;199:108012. doi:10.1016/j.csda.2024.108012. [Google Scholar] [CrossRef]

6. Wang X, Chen H, Cai W, Shen D, Huang H. Regularized modal regression with applications in cognitive impairment prediction. Adv Neural Inf Process Syst. 2017;30:1448–58. [Google Scholar] [PubMed]

7. Menezes AF, Mazucheli J, Chakraborty S. A collection of parametric modal regression models for bounded data. J Biopharm Stat. 2021;31(4):490–506. doi:10.1080/10543406.2021.1918141. [Google Scholar] [PubMed] [CrossRef]

8. Hassan AS, Metwally DS, Elgarhy M, Semary HE, Faal A, Mohamed RE. Sine power unit inverse lindley model: Bayesian analysis and practical application. Eng Rep. 2025;7(6):e70242. [Google Scholar]

9. Muhammad M, Abba B, Xiao J, Alsadat N, Jamal F, Elgarhy M. A new three-parameter flexible unit distribution and its quantile regression model. IEEE Access. 2024;12(149):156235–51. doi:10.1109/access.2024.3485219. [Google Scholar] [CrossRef]

10. Elgarhy M, Abdalla GSS, Hassan AS, Almetwally EM. Bayesian and non-bayesian analysis of the novel unit inverse exponentiated Lomax distribution using progressive censoring schemes with optimal scheme and data application. Comput J Math Stat Sci. 2025;5(1):78–108. doi:10.21608/cjmss.2025.374277.1151. [Google Scholar] [CrossRef]

11. Ghosh I. Bradford distribution and its application in modeling medical data: a suitable alternative to distributions defined on the unit interval. J Appl Stat. 2026;53(3):520–36. doi:10.1080/02664763.2025.2520342. [Google Scholar] [PubMed] [CrossRef]

12. Whittaker ET, Watson GN. A course of modern analysis. Cambridge mathematical library. 4th ed. Cambridge, UK: Cambridge University Press; 2020. [Google Scholar]

13. Dunn PK, Smyth GK. Randomized quantile residuals. J Comput Graph Stat. 1996;5(3):236–44. doi:10.1080/10618600.1996.10474708. [Google Scholar] [CrossRef]

14. Nasiru S, Chesneau C, Ocloo SK. The log-cosine-power unit distribution: a new unit distribution for proportion data analysis. Decis Anal J. 2024;10:100397. [Google Scholar]

15. Krishna A, Maya R, Chesneau C, Irshad MR. The unit Teissier distribution and its applications. Math Comput Appl. 2022;27(1):12. doi:10.3390/mca27010012. [Google Scholar] [CrossRef]

16. Korkmaz MC, Chesneau C, Korkmaz ZS. Transmuted unit Rayleigh quantile regression model: alternative to beta and Kumaraswamy quantile regression models. Univ Politeh Buchar Sci Bull Ser Appl Math Phys. 2021;83(3):149–58. [Google Scholar]

17. Topp CW, Leone FC. A family of J-shaped frequency functions. J Am Stat Assoc. 1955;50(269):209–19. doi:10.1080/01621459.1955.10501259. [Google Scholar] [CrossRef]

18. Bakouch HS, Hussain T, Tosic M, Stojanovic VS, Qarmalah N. Unit exponential probability distribution: characterization and applications in environmental and engineering data modeling. Mathematics. 2023;11(19):4207. [Google Scholar]

19. Ribeiro TF, Pena-Ramirez FA, Guerra RR, Cordeiro GM. Another unit Burr XII quantile regression model based on the different reparameterization applied to dropout in Brazilian undergraduate courses. PLoS One. 2022;17(11):e0276695. doi:10.1371/journal.pone.0276695. [Google Scholar] [PubMed] [CrossRef]

20. Kumaraswamy P. A generalized probability density function for double-bounded random processes. J Hydrol. 1980;46(1–2):79–88. doi:10.1016/0022-1694(80)90036-0. [Google Scholar] [CrossRef]