The Continuous Bernoulli Distribution: Mathematical Characterization, Fractile Regression, Computational Simulations, and Applications

Abstract

The continuous Bernoulli distribution is defined on the unit interval and has a unique property related to fractiles. A fractile is a position on a probability density function where the corresponding surface is a fixed proportion. This article presents the derivation of properties of the continuous Bernoulli distribution and formulates a fractile or quantile regression model for a unit response using the exponentiated continuous Bernoulli distribution. Monte Carlo simulation studies evaluate the performance of point and interval estimators for both the continuous Bernoulli distribution and the fractile regression model. Real-world datasets from science and education are analyzed to illustrate the modeling abilities of the continuous Bernoulli distribution and the exponentiated continuous Bernoulli quantile regression model.

1. Introduction

The interest in new probability discrete and continuous distributions is increasing; see, for example, [1]. Recently, this interest has been present in distributions defined on the interval [0, 1], which has raised great attention from researchers. These distributions have many applications in several fields of science; see, for example, [2,3,4,5,6,7]. Often, novel unit models are transported to unit intervals by transformations of well-known models [2], with studies on unit modeling increasing day by day [2,3,4,5,6,7].

Regression is a statistical approach that is widely employed [8] to explain the structure between a dependent (response) variable and independent (explanatory) variables or covariates. Such a structure can be linear or nonlinear, associating the response mean with the independent variables. When the dependent variable is stated on the interval [0, 1], a regression structure employing the beta distribution [9] has been the most preferred method to explain the unit-dependent variable with the independent variables based on the conditional mean. Note that the beta model is formulated in relation to its mean. Then, its unit regression based on the mean was established in [9].

With similar ideas,  alternative unit mean response regressions were introduced, among them are the beta rectangular [10], Birnbaum–Saunders [11], log-Bilal [12], log-Lindley [1], log-weighted exponential [13], and unit Lindley [14] models.

The dependent variables can be affected by situations related to a atypical data or asymmetric distributions. When these situations exist, the mean is also affected strongly, and the inference from the model may give rise to possible inaccurate interpretations.

A fractile is a point on a probability density function (PDF) where the surface upon the function is a fixed proportion. Due to the inaccuracy mentioned, the quantile regression (QR) models are a robust alternative to mean regression models. The QR model was introduced in [15] and relates the response’s conditional quantiles to independent variables’ given values instead of explaining their conditional mean. One of the key benefits of using QR models over standard mean regressions is its lack of imposed distributional assumptions on the error term. This allows for greater flexibility in modeling different types of data [16].

To implement a QR model utilizing a probability distribution, it is necessary to parameterize it using its quantile function (QF). Fortunately, this parameterization can be utilized for any probability distribution with a closed-form QF, regardless of whether its mean has a closed form. This provides a high level of adaptability and enables us to model quantiles based on various distributions, as demonstrated in the following examples: the Kumaraswamy (KW) [17], L-logistic [18], log-extended exponential-geometric (LEEG) [19], log-symmetric [20], power Johnson SB, [21], unit-Birnbaum–Saunders [22], unit-Burr-XII, [2,23], unit-Chen [24], unit-Weibull [25], and Weibull–Marshall–Olkin [26] QR models. These types of models and structures have been largely applied to describe COVID-19 mortality [23,27].

The continuous Bernoulli (CB) model, as described in [4], is a probability distribution that exists on the interval [0, 1] and has a single parameter related to the shape of the distribution. It belongs to the exponential family of distributions and offers a generalization of the uniform distribution defined on the unit interval. In their study, the authors of [4] proposed using the CB distribution to investigate the impact of Bernoulli variational autoencoders on intensity data that fall within the [0, 1] range. The results of their research indicated that the CB distribution leads to significant improvements in various metrics and datasets, including the generation of sharper image samples, and suggested a wider range of potential performance-optimized variational autoencoders. Despite the promising findings of the CB distribution, as far as we know, its associated mathematical and statistical properties as well as its associated QR model have yet to be explored. Thus, the primary objectives of this article are as follows:

(i)

To derive some mathematical and statistical properties that have not been studied of the CB distribution stated on the interval [0, 1], including estimation and inference.

(ii)

To formulate a QR model for a unit response employing the exponentiated CB (ECB) distribution. In this way, the applications of the CB distribution are extended. The formulated regression model is a robust alternative to the beta regression and other models proposed in the literature. Another motivation for the proposed QR model comes from a real-world dataset.

The present article is organized into seven main sections. Section 2 is dedicated to deriving mathematical properties of the continuous Bernoulli (CB) distribution, while Section 3 focuses on the various point estimation procedures for its parameters. A novel QR model is introduced in Section 4. Monte Carlo simulations are used to evaluate the performance of point and interval estimators for both the CB distribution and the associated QR model in Section 5. The results are illustrated utilizing a real-world dataset in Section 6. Finally, Section 7 sketches some concluding remarks.

2. The Continuous Bernoulli Distribution

This section provides some known and new functions of the CB distribution.

2.1. Cumulative Distribution and Probability Density Functions

In the fields of artificial intelligence, in the context of modeling pixel intensities, the CB distribution is a commonly encountered concept [28]. The PDF and cumulative distribution function (CDF) of the CB distribution are formulated as

$$F(x,\beta )=\left\{\begin{array}{cc}x,\hfill & \beta =1/2;\hfill \\ \left(1-\beta \right)\left({\left({\displaystyle \frac{\beta }{1-\beta }}\right)}^x-1\right){\displaystyle \frac{1}{\left(2\beta -1\right)}},\hfill & \beta \ne 1/2;\hfill \end{array}\right.$$

(1)

and

$$f(x,\beta )=\left\{\begin{array}{cc}1,\hfill & \beta =1/2;\hfill \\ \left(1-\beta \right){\left({\displaystyle \frac{\beta }{1-\beta }}\right)}^x{\displaystyle \frac{1}{\left(2\beta -1\right)}}log\left(\beta /(1-\beta )\right),\hfill & \beta \ne 1/2;\hfill \end{array}\right.$$

respectively, with $x\in (0,1)$, and $\beta \in \left(0,1\right)$ being a parameter related to the shape of the distribution. Alternatively, the PDF of the CB distribution can be written as

$$f\left(x,\beta \right)=\left\{\begin{array}{cc}1,\hfill & \beta =1/2;\hfill \\ {\displaystyle \frac{2\left(1-\beta \right){tanh}^{-1}\left(1-2\beta \right)}{\left(1-2\beta \right)}}{\left(\beta /(1-\beta )\right)}^x,\hfill & \beta \ne 1/2;\hfill \end{array}\right.$$

where ${tanh}^{-1}(w)=(1/2)log((1+w)/(1-w))$, for $w\in (-1,1)$, is the inverse hyperbolic tangent function. We denote the CB distribution of parameter $\beta$ as CB$(\beta )$. Note that the CB distribution corresponds to the uniform distribution on the interval [0, 1], that is, on (0, 1), when $\beta =1/2$. Observe that there is an interesting recursive relation between the CB PDF and its derivatives. For $\beta \ne 1/2$, this relation is stated as

$${\displaystyle \frac{{\mathrm{d}}^{m}f\left(x,\beta \right)}{\mathrm{d}{x}^m}}=f(x,\beta ){log}^{m+1}\left({\displaystyle \frac{\beta }{1-\beta }}\right),m\in \{1,2,\cdots \}.$$

Hence, it can be concluded that the PDF with a concave CDF is decreasing for $\beta <1/2$, and the PDF with a convex CDF is increasing for $\beta >1/2$. Note that the term ${\sum }_{i=1}^n{X}_i$ is a sufficient statistic for estimating $\beta$ based on a sample of size n, $X_1,\cdots ,X_n$ namely, collected from the CB distribution.

2.2. Inverse Cumulative Distribution Function

From (1), we can obtain the CB inverse CDF, which is given by

$$F^{-1}(x,\beta )=\left\{\begin{array}{cc}x,\hfill & \beta =1/2;\hfill \\ {\displaystyle \frac{log\left(1+\left({\displaystyle \frac{2\beta -1}{1-\beta }}\right)x\right)}{log\left({\displaystyle \frac{\beta }{1-\beta }}\right)}},\hfill & \beta \ne 1/2.\hfill \end{array}\right.$$

(2)

In Figure 1, we can see the plot of (2) for the values $\beta =0.25$, $\beta =0.5$, and $\beta =0.75$. To simulate a random variable (RV) with CB distribution, we propose Algorithm 1. In simulation studies related to the CB distribution, generating data from this distribution is important. Then, we can use Algorithm 1 for this purpose, which takes advantage of (2) to generate a sample of size n.

Algorithm 1: Approach that simulates the CB distribution using the inverse-transform method

Computational experiments were conducted utilizing Algorithm 1.

Table 1. Simulation results with $\beta =0.25$.

n	Runtime (in Seconds)	$\overline{\mathit{x}}$	${\mathit{s}}^2$
25	$0.0116$	$0.4106$	$0.0780877$
50	$0.0207$	$0.4110$	$0.0788680$
100	$0.0352$	$0.4101$	$0.0783473$
200	$0.0513$	$0.4098$	$0.0785829$
500	$0.1273$	$0.4104$	$0.0784173$
1000	$0.1951$	$0.4102$	$0.0784991$

displays the results generated with $\beta =0.25$, while

Table 2. Simulation results with $\beta =0.75$.

n	Runtime (in Seconds)	$\overline{\mathit{x}}$	${\mathit{s}}^2$
25	$0.0197$	$0.5895$	$0.0784745$
50	$0.0224$	$0.5897$	$0.0780702$
100	$0.0315$	$0.5885$	$0.0788916$
200	$0.0459$	$0.5887$	$0.0786282$
500	$0.0955$	$0.5894$	$0.0788145$
1000	$0.2136$	$0.5897$	$0.0784867$

provides the results generated with $\beta =0.75$, using $n\in \{25,50,100,200,500,1000\}$ and performing 1000 Monte Carlo simulations for each n.

The simulation results demonstrate the experimental convergence of Algorithm 1. The theoretical mean and variance values for the CB distribution with $\beta =0.25$ are $0.4102$ and $0.0785354$, respectively, while the theoretical mean and variance values with $\beta =0.75$ are $0.5898$ and $0.0785354$, respectively.

2.3. Hazard Rate Function

The hazard rate function (HRF) of the CB distribution is stated as

$$h(x,\beta )=\frac{{\left(1-\beta \right)}^{1-x}{\beta }^x}{\beta -{\left(1-\beta \right)}^{1-x}{\beta }^x}log\left(\beta /(1-\beta )\right),x\in (0,1),\beta \ne 1/2.$$

The HRF is an increasing function in x because its derivative expressed as

$$\frac{\mathrm{d}h\left(x,\beta \right)}{\mathrm{d}x}=\frac{\beta \left(1-\beta \right)}{{\left(\beta -\left(1-\beta \right){\left(\beta /(1-\beta )\right)}^x\right)}^2}{\left(\beta /(1-\beta )\right)}^x{log}^2\left(\beta /(1-\beta )\right)$$

is positive for all $\beta$. Furthermore, note that ${lim}_{x\to 1}h(x,\beta )=\infty$ and

$$\underset{x\to 0}{lim}h\left(x,\beta \right)=\left\{\begin{array}{cc}1,\hfill & \beta =1/2;\hfill \\ \frac{\left(1-\beta \right)}{2\beta -1}log\left(\beta /(1-\beta )\right),\hfill & \beta \ne 1/2.\hfill \end{array}\right.$$

Figure 2 displays the CB PDF, CDF, and HRF. From this figure, note the flexible shapes these functions have for the different parameter values, that is, for $\beta$.

2.4. Moments

Using the definition of the moment-generating function, for the case of the CB distribution, we obtain that

$$M_X\left(t\right)=\mathrm{E}\left({\mathrm{e}}^{tX}\right)=\frac{\left(\beta \left(1+{\mathrm{e}}^t\right)-1\right)log\left(\beta /(1-\beta )\right)}{\left(2\beta -1\right)\left(log\left(\beta /(1-\beta )\right)+t\right)},\beta \ne 1/2,t\in \Re .$$

Note that $M_X(t)=({\mathrm{e}}^t-1)/t$ for $\beta =1/2$. Thus, utilizing the relation given by $\mathrm{E}(X^r)={\mathrm{d}}^r{M}_X(t)/\mathrm{d}{t}^r{|}_{t=0}$, for $r\in \{1,2,\cdots \}$, the first four noncentral moments are, respectively, for $\beta \ne 1/2$, stated as

$$\begin{array}{ccc}\hfill m_1=\mathrm{E}(X)& =& \frac{\beta log\left(\beta /(1-\beta )\right)+1-2\beta }{\left(2\beta -1\right)log\left(\beta /(1-\beta )\right)},\hfill \\ \hfill m_2=\mathrm{E}(X^2)& =& \frac{\beta {log}^2\beta -2\beta log\beta log\left(1-\beta \right)+\beta {log}^2\left(1-\beta \right)-2\beta log\beta +2\beta log\left(1-\beta \right)+4\beta -2}{\left(2\beta -1\right){log}^2\left(\beta /(1-\beta )\right)},\hfill \\ \hfill m_3=\mathrm{E}(X^3)& =& \frac{-\beta {log}^3\left(1-\beta \right)+3\beta \left(log\beta -1\right){log}^2\left(1-\beta \right)-3\beta \left({log}^2\beta +2-2log\beta \right)log\left(1-\beta \right)}{\left(2\beta -1\right){log}^3\left(\beta /(1-\beta )\right)}\hfill \\ \hfill & & +\frac{\beta {log}^3\beta -12\beta +6\beta log\beta -3\beta {log}^2\beta +6}{\left(2\beta -1\right){log}^3\left(\beta /(1-\beta )\right)},\hfill \\ \hfill m_4=\mathrm{E}(X^4)& =& \frac{\beta {log}^4\beta -4\beta {log}^3\beta log\left(1-\beta \right)+6\beta {log}^2\beta {log}^2\left(1-\beta \right)-4\beta {log}^3\left(1-\beta \right)log\beta }{\left(2\beta -1\right){log}^4\left(\beta /(1-\beta )\right)}\hfill \\ \hfill & & +\frac{\beta {log}^4\left(1-\beta \right)-4\beta {log}^3\beta +12\beta {log}^2\beta log\left(1-\beta \right)-12\beta {log}^2\left(1-\beta \right)log\beta +4\beta {log}^3\left(1-\beta \right)+12\beta {log}^2\beta }{\left(2\beta -1\right){log}^4\left(\beta /(1-\beta )\right)}\hfill \\ \hfill & & +\frac{-24\beta log\beta log\left(1-\beta \right)+12\beta {log}^2\left(1-\beta \right)-24\beta log\beta +24\beta log\left(1-\beta \right)+48\beta -24}{\left(2\beta -1\right){log}^4\left(\beta /(1-\beta )\right)}.\hfill \end{array}$$

Figure 3 shows some shapes of the CB moments. Note that, as $\beta$ increases, the mean increases, whereas the skewness decreases. According to the values of $\beta$, the variance increases for $\beta <1/2$ and decreases for $\beta >1/2$, but a reverse situation is detected for the kurtosis.

2.5. Quantile Function

The QF, defined by $F^{-1}(q,\beta )=x_q(\beta )$, for $q\in (0,1)$, of the CB distribution can be obtained in a closed form. Let X be a CB distributed RV with CDF as defined in (1). Then, the CB QF is given by

$$x_q\left(\beta \right)=\frac{log\left(1+\frac{\left(2\beta -1\right)q}{1-\beta }\right)}{log\left(\beta /(1-\beta )\right)},\beta \ne 1/2,q\in (0,1).$$

(3)

Note that the CB QF is equal to q for $\beta =1/2$, whereas the corresponding median is derived as

$$x_{0.5}(\beta )={\left(log\left(\beta /(1-\beta )\right)\right)}^{-1}log\left(\frac{0.5}{1-\beta }\right),$$

for $\beta \ne 1/2$. Notice that if the RV Q related to q in (3) has a uniform distribution on (0, 1), which is denoted as U(0, 1), then the RV associated with $x_Q(\beta )$ has the CB distribution.

2.6. Mean Residual Life Function

The mean residual life (MRL), also known as life expectancy, is a critical parameter in the field of reliability that characterizes the assumed distribution. Unlike the hazard rate function (HRF), which only considers the risk of immediate failure, the MRL summarizes the entire residual life distribution. The MRL function of a positive random variable at time x is defined as

$$\mathrm{MRL}(x)=\frac{1}{1-F(x)}{\int }_x^{\infty }\left(1-F\left(y\right)\right)\mathrm{d}y,$$

where F is the lifetime model CDF. Therefore, the MRL function of the CB distribution is stated as

$$\mathrm{MRL}(x)=\frac{\beta \left(1-x\right)log\left(\beta /(1-\beta )\right)-\beta \left(1+{\left(\beta /(1-\beta )\right)}^x\right)+{\left(\beta /(1-\beta )\right)}^x}{\left(\beta -\left(1-\beta \right){\left(\beta /(1-\beta )\right)}^x\right)log\left(\beta /(1-\beta )\right)},\beta \ne 1/2,$$

while MRL(0) is obtained as

$$\mathrm{MRL}(0)=\frac{\beta log\left(\beta /(1-\beta )\right)+1-2\beta }{\left(2\beta -1\right)log\left(\beta /(1-\beta )\right)},\beta \ne 1/2,$$

which is equal to the corresponding expected value, and then MRL(1) = 0. Observe that $\mathrm{MRL}(x)=(1-x)/2$, for $\beta =1/2$.

We can comment about the shapes of the MRL function via its HRF. The following results hold for both discrete and continuous lifetimes [29]:

• If the HRF is (strictly) increasing, then the MRL is (strictly) decreasing.
• If the HRF is (strictly) decreasing, then the MRL is (strictly) increasing.
• The HRF is a constant function (exponential or geometric distribution) if and only if the MRL is a constant.

We can conclude that, since the HRF function of the CB distribution is strictly increasing, its MRL function is strictly decreasing. For detail about the MRL, see [30].

2.7. Incomplete Moments and Lorenz Curve

For $\beta \ne 1/2$ and $r\in \{1,2,\cdots \}$, the rth incomplete moment of the CB distribution is given by

$$\begin{array}{ccc}\hfill I(x,r)& =& {\int }_0^x{z}^{r}f(z,\beta )\mathrm{d}z\hfill \\ & =& \frac{x^r\left(1-\beta \right)\left(r\mathsf{\Gamma }\left(r,-xlog\left(\beta /(1-\beta )\right)\right){\left(-xlog\left(\beta /(1-\beta )\right)\right)}^{-r}-\mathsf{\Gamma }\left(r+1\right){\left(-xlog\left(\beta /(1-\beta )\right)\right)}^{-r}+{\left(\beta /(1-\beta )\right)}^x\right)}{2\beta -1},\hfill \end{array}$$

for $z<x$, where $\mathsf{\Gamma }(·)$ and $\mathsf{\Gamma }(·,·)$ are the gamma function and upper incomplete gamma function, respectively. For $\beta =1/2$, the CB rth incomplete moment is equal to $x^{r+1}/(r+1)$. Note that the equations ${lim}_{x\to \infty }I(x,r)=\mathrm{E}(X^r)$ and $F(y,\beta )=I(x,r=0)$ hold.

As a result of the incomplete moments, the Lorenz curve of the CB distribution for $\beta \ne 1/2$ is expressed by

$$\begin{array}{ccc}\hfill L(p)& =& \frac{1}{\mathrm{E}(X)}{\int }_0^{x_q(\beta )}xf(x,\beta )\mathrm{d}x\hfill \\ & =& \frac{I\left(x_q\left(\beta \right),1\right)}{\mathrm{E}(X)}=\frac{\left(\beta \left(2p-1\right)+1-p\right)log\left(\frac{1+\beta \left(2p-1\right)-p}{1-\beta }\right)-2p\beta +p}{1-2\beta +\beta log\left(\beta /(1-\beta )\right)}.\hfill \end{array}$$

Then, $L(p)=p^2$ is obtained for $\beta =1/2$.

2.8. Stress Strength Reliability

The system with applied stress works smoothly if its strength is higher than the applied stress. The stress–strength parameter of a life distribution is an important indicator of reliability. This parameter is the working probability of the system and a reliability measure, which is defined as

$$R=P\left(X<Y\right)={\int }_0^1{f}_X\left(x\right)F_Y\left(x\right)\mathrm{d}x,$$

where $f_X$ is the PDF of the stress (X) and $F_Y$ is the CDF of the strength (Y). Now, consider X and Y as two independent CB distributed RVs with parameters ${\beta }_1$ and ${\beta }_2$, respectively. Then, the CB stress–strength parameter is given by

$$R=\frac{\left(2{\beta }_1\left({\beta }_2-1\right)-{\beta }_2+1\right)log\left(\frac{{\beta }_2}{1-{\beta }_2}\right)+\left(2{\beta }_2-1\right){\beta }_{1}log\left(\frac{{\beta }_1}{1-{\beta }_1}\right)}{\left(2{\beta }_2-1\right)\left(2{\beta }_1-1\right)log\left(\frac{{\beta }_1}{(1-{\beta }_1)}\frac{{\beta }_2}{(1-{\beta }_2)}\right)}.$$

Note that, when ${\beta }_1={\beta }_2\ne 1/2$, then $R=1/2$.

3. Point Estimation for the Parameter of the CB Distribution

This section provides some point estimation methods, such as maximum likelihood (ML), moment, percentile (PE), least square (LS), weighted square (WS), Anderson–Darling (AD), and Cramér–von Mises (CM) to determine the model parameter $\beta$. Note that, for the estimation processes of the above methods, we focus on $\beta \ne 1/2$.

3.1. Maximum Likelihood Method

Consider $X_1,\cdots ,X_n$ as a sample of size n collected from the CB distribution, with $x_1,\cdots ,x_n$ as its observed values. Thus, the associated log-likelihood function is established as

$$\begin{array}{cc}\hfill \ell \left(\beta \right)& =nlog\left(1-\beta \right)-nlog\left(2\beta -1\right)+nlog\left(log\left(\beta /(1-\beta )\right)\right)+log\left(\beta /(1-\beta )\right)\sum _{i=1}^n{x}_i.\hfill \end{array}$$

(4)

Taking the first derivative of the expression stated in (4) with respect to $\beta$, we obtain

$$\frac{\mathrm{d}\ell (\beta )}{\mathrm{d}\beta }=-\frac{n}{log\left(1-\beta \right)}-\frac{2n}{log\left(2\beta -1\right)}+\frac{n}{\beta \left(1-\beta \right)log\left(\beta /(1-\beta )\right)}+\frac{1}{\beta \left(1-\beta \right)}\sum _{i=1}^n{x}_i.$$

(5)

Since the equation $\mathrm{d}\ell (\beta )/\mathrm{d}\beta =0$ obtained from (5) is a nonlinear function of $\beta$, the ML estimator $\widehat{\beta }$ must be computed using numerical methods. For any smooth PDF, it is widely recognized that under regularity conditions, $\widehat{\beta }$ has a normal distribution asymptotically with mean $\beta$ and variance ${\mathcal{I}}_{\beta }^{-1}$, with ${\mathcal{I}}_{\beta }=-\mathrm{E}({\mathrm{d}}^2\ell (\beta )/\mathrm{d}{\beta }^2)$ being the Fisher information.

After performing several computations, the inverse of the Fisher information can be obtained as

$${\mathcal{I}}_{\beta }^{-1}=\frac{{\left(\beta \left(1-\beta \right)\left(2\beta -1\right)log\left(\beta /(1-\beta )\right)\right)}^2}{n\left(1-\beta \left(1-\beta \right)\left(4+{log}^2\left(\beta /(1-\beta )\right)\right)\right)}.$$

(6)

Based on the previous results, the asymptotic $100(1-\theta )\%$ confidence interval of the parameter $\beta$ is given as

$$\widehat{\beta }-z_{\theta /2}\frac{\beta \left(1-\beta \right)\left(2\beta -1\right)log\left(\beta /(1-\beta )\right)}{\sqrt{n\left(1-\beta \left(1-\beta \right)\left(4+{log}^2\left(\beta /(1-\beta )\right)\right)\right)}}\le \beta \le \widehat{\beta }+z_{\theta /2}\frac{\beta \left(1-\beta \right)\left(2\beta -1\right)log\left(\beta /(1-\beta )\right)}{\sqrt{n\left(1-\beta \left(1-\beta \right)\left(4+{log}^2\left(\beta /(1-\beta )\right)\right)\right)}},$$

where $z_{\theta /2}$ refers to the $(\theta /2)\times 100$th upper percentile of the standard normal distribution.

3.2. Moment Estimation

The parameter estimation for the CB distribution can be performed using the moment method, as its expected value can be expressed analytically. This method involves equating the first moment of the sample to the expected value of the CB distribution. By solving this equation, the moment estimate of the parameter $\beta$ can be obtained as

$$\overline{x}-\frac{\beta log\left(\beta /(1-\beta )\right)+1-2\beta }{\left(2\beta -1\right)log\left(\beta /(1-\beta )\right)}=0,$$

(7)

where $\overline{x}$ is the sample mean. Note that the equation stated in (7) must be solved employing a numerical method.

3.3. Percentile Estimation

The percentile method to estimate $\beta$ in the CB distribution was originally proposed in [31]. Since the QF of the CB distribution has a closed form, this method can be applied straightforwardly. Consider $X_{(1)},\cdots ,X_{(n)}$ as the order statistics of a sample of size n from the CB distribution and $x_{(1)},\cdots ,x_{(n)}$ as their observed values. An estimate of the CDF can be obtained as $\widehat{F}(x_{(i)})=i/(n+1)$, for $i\in \{1,\cdots ,n\}$. Then, an estimate of $\beta$ can be reached by minimizing a certain function with respect to $\beta$. The expression for this function is given by

$$\begin{array}{cc}\hfill \mathrm{PE}(\beta )& =\sum _{i=1}^n{\left(x_{\left(i\right)}-x_{q_i}\left(\beta \right)\right)}^2,\hfill \end{array}$$

(8)

and then the PE estimate of $\beta$ is attained, where $x_q$ is defined in (3). Therefore, the PE estimate of $\beta$ is the solution on $\beta$ evaluated at $\widehat{\beta }$ of an equation formulated as

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}\mathrm{PE}(\beta )}{\mathrm{d}\beta }& =& \frac{2}{\beta \left(1-\beta \right){log}^2\left(\beta /(1-\beta )\right)}{\displaystyle \sum _{i=1}^n}\left(x_{\left(i\right)}-\frac{log\left(1+\frac{i\left(2\beta -1\right)}{\left(n+1\right)\left(1-\beta \right)}\right)}{log\left(\beta /(1-\beta )\right)}\right)\hfill \\ & & \times \frac{\left(\beta \left((2i-1-n)\right)+n-i+1\right)log\left(1+\frac{i\left(2\beta -1\right)}{\left(n+1\right)\left(1-\beta \right)}\right)-i\beta log\left(\beta /(1-\beta )\right)}{\left(n+1\right)\left(1-\beta \right)+i\left(2\beta -1\right)}=0.\hfill \end{array}$$

3.4. Other Methods of Estimation

Next, the LS, WS, AD, and CM methods of estimation for $\beta$ are provided. According to these methods, the functions to be used for estimating $\beta$ are expressed as

$$\begin{array}{ccc}\hfill \mathrm{LS}\left(\beta \right)& =& {\displaystyle \sum _{i=1}^n}{\left(F(x_{(i)},\beta )-\frac{i}{n+1}\right)}^2,\hfill \\ \hfill \mathrm{WS}\left(\beta \right)& =& {\displaystyle \sum _{i=1}^n}\frac{(n+2){(n+1)}^2}{(n-i+1)i}{\left(F(x_{(i)},\beta )-\frac{i}{n+1}\right)}^2,\hfill \\ \hfill \mathrm{LS}\left(\beta \right)& =& -n-{\displaystyle \sum _{i=1}^n}\frac{(2i-1)}{n}\left(log\left(1-F(x_{(n+1-i)},\beta )\right)+log(F(x_{(i)},\beta ))\right),\hfill \\ \hfill \mathrm{CM}\left(\beta \right)& =& \frac{1}{12n}+{\displaystyle \sum _{i=1}^n}{\left(F(x_{(i)},\beta )-\frac{2i-1}{2n}\right)}^2,\hfill \end{array}$$

In the CB distribution, a sample of size n is drawn from the population, with $x_{(1)},\cdots ,x_{(n)}$ representing the ordered observed values of this sample. Therefore, the LS, WS, AD, and CM estimates of $\beta$ are given by

$$\begin{array}{ccc}\hfill {\widehat{\beta }}_{\mathrm{LS}}& =& \underset{\beta }{argmin}\left\{\mathrm{LS}\left(\beta \right)\right\},\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\widehat{\beta }}_{\mathrm{WS}}& =& \underset{\beta }{argmin}\left\{\mathrm{WS}\left(\beta \right)\right\},\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill {\widehat{\beta }}_{\mathrm{AD}}& =& \underset{\beta }{argmin}\left\{\mathrm{AD}\left(\beta \right)\right\},\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\widehat{\beta }}_{\mathrm{CM}}& =& \underset{\beta }{argmin}\left\{\mathrm{CM}\left(\beta \right)\right\},\hfill \end{array}$$

(12)

respectively. It is important to emphasize that the equations presented in (4), (7), (8), (9), (10), (11) and (12) do not have explicit solutions and must be solved through numerical methods. However, these equations can be optimized directly using software such as R or Matlab.

4. Associated QR Model

This section states a QR model employing the exponentiated CB distribution.

4.1. The Exponentiated CB Distribution

Although the CB distribution has a closed-form mean function and QF, its PDF and CDF cannot be parameterized using these functions due to the dependence of the parameter $\beta$. Therefore, we propose an alternative QR model that is based on the ECB distribution, which is a generalized form of the CB distribution. The ECB CDF and PDF are expressed as

$$\begin{array}{c}\hfill G\left(y,\alpha ,\beta \right)=\left\{\begin{array}{cc}{y}^{\alpha },\hfill & \beta =1/2;\hfill \\ {\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^y-1\right)}{\left(2\beta -1\right)}\right)}^{\alpha },\hfill & \beta \ne 1/2;\hfill \end{array}\right.\end{array}$$

and

$$\begin{array}{c}\hfill g\left(y,\alpha ,\beta \right)=\left\{\begin{array}{cc}\alpha y^{\alpha -1},\hfill & \beta =1/2;\hfill \\ \frac{\alpha {\left(\beta /(1-\beta )\right)}^{y}log\left(\beta /(1-\beta )\right)}{{\left(\beta /(1-\beta )\right)}^y-1}{\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^y-1\right)}{\left(2\beta -1\right)}\right)}^{\alpha },\hfill & \beta \ne 1/2;\hfill \end{array}\right.\end{array}$$

respectively, where $y\in (0,1)$, $\alpha >0$ is a parameter related to the shape of the distribution and $\beta \in (0,1)$. The CDF of the ECB distribution is $G(y,\alpha ,\beta )={(F(y,\beta ))}^{\alpha }$, and it can be named a Lehmann-type I CB distribution. For $\beta =1/2$, the ECB distribution is the power function distribution.

The $\tau \times 100$th quantile of the ECB distribution is presented by

$$Q(\tau ,\alpha ,\beta )=\left\{\begin{array}{cc}{\tau }^{1/\alpha },\hfill & \beta =1/2;\hfill \\ {\left(log\left(\beta /(1-\beta )\right)\right)}^{-1}log\left(1+\frac{\left(2\beta -1\right){\tau }^{1/\alpha }}{1-\beta }\right),\hfill & \beta \ne 1/2;\hfill \end{array}\right.$$

where $\tau \in (0,1)$.

By using the QF of the ECB distribution, we can now reparameterize it in terms of its quantiles. Observe that, as $\beta =1/2$, the ECB distribution can be reparameterized as a particular case of the Kumaraswamy distribution based on the quantile form, as presented in [17]. Here, we consider $\beta \ne 1/2$ and provide the PDF form of the ECB distribution based on its QF.

Consider $\mu =Q(\tau ,\alpha ,\beta )$ and $\alpha =log(\tau {\left(log\left(\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{\mu }-1\right)/\left(2\beta -1\right)\right)\right)}^{-1})$. Hence, the reparameterized ECB CDF and PDF are introduced by

$$G\left(y,\beta ,\mu ,\tau \right)={\tau }^{log\left(\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^y-1\right)/\left(2\beta -1\right)\right){\left(log\left(\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{\mu }-1\right)/\left(2\beta -1\right)\right)\right)}^{-1}}$$

(13)

and

$$\begin{array}{c}\hfill g\left(y,\beta ,\mu ,\tau \right)=\frac{log\tau {\left(\beta /(1-\beta )\right)}^{y}log\left(\beta /(1-\beta )\right){\tau }^{log\left(\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^y-1\right)/\left(2\beta -1\right)\right){\left(log\left(\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{\mu }-1\right)/\left(2\beta -1\right)\right)\right)}^{-1}}}{log\left(\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{\mu }-1\right)/\left(2\beta -1\right)\right)\left({\left(\beta /(1-\beta )\right)}^y-1\right)}.\end{array}$$

(14)

The reparameterized PDF stated in (14) is defined by the CDF presented in (13), where $\beta \in (0,1)$ is a parameter related to the shape of the distribution, $\mu \in (0,1)$ is the quantile parameter, and $\tau \in (0,1)$ is a known quantile level. We denote the RV associated with the CDF given in (13) as $Y\sim \mathrm{ECB}(\beta ,\mu ,\tau )$. Figure 4 illustrates the diverse shapes of the reparameterized PDF, which exhibit characteristics such as U-shape, unimodality, increasing, decreasing, and constancy.

4.2. ECB Quantile Regression Model

Once the reparameterized ECB distribution is defined, a QR model can be introduced. Suppose we have n observations $y_1,\cdots ,y_n$ from the reparameterized distribution with its PDF stated in (14). Then, the ECB QR model can be expressed as

$$g({\mu }_i)={\mathbf{x}}_i^{\top }\mathbf{\beta }.$$

(15)

Here, ${\mathbf{x}}_i={\left(1,x_{i1},\cdots ,x_{ip}\right)}^{\top }$ is the vector of independent variable values, $\mathbf{\beta }={\left({\beta }_0,{\beta }_1,\cdots ,{\beta }_p\right)}^{\top }$ is the vector of unknown regression coefficients, and $g:\left(0,1\right)\to \Re$ is a strictly monotonic and twice differentiable link function used to relate independent variables with the conditional quantile of the dependent variable. Common link functions include the logit, probit, log-log, and Cauchy structures. Thus, ${\mu }_i$ is obtained by inverting the model given in (15) such that ${\mu }_i=g^{-1}({\mathbf{x}}_i^{\top }\mathbf{\beta })$, for $i\in \{1,\cdots ,n\}$. In this proposed regression model, we use the logit function for linking, stated as $g({\mu }_i)=log({\mu }_i/(1-{\mu }_i))$, for $i\in \{1,\cdots n\}$. If we take $\tau =0.5$, the link function connects the median response to the independent variables in the ECB QR model.

4.3. ML Estimation of the Regression Coefficients

Assume that we have n independent RVs $Y_1,\cdots ,Y_n$ following the ECB distribution, with $\beta$ being a parameter related to the shape of the distribution, ${\mu }_i$ a quantile parameter, and $\tau$ a known quantile level. The quantile parameter is obtained using the logistic function formulated as

$${\mu }_i=\frac{{\mathrm{e}}^{{\mathbf{x}}_i^{\top }\mathbf{\beta }}}{1+{\mathrm{e}}^{{\mathbf{x}}_i^{\top }\mathbf{\beta }}},$$

(16)

where ${\mathbf{x}}_i$ is the vector of independent variable values for the ith observation, and $\mathbf{\beta }$ is the vector of unknown regression coefficients. The parameter vector $\mathbf{\vartheta }=(\beta ,{\mathbf{\beta }}^{\top })$ has $p+2$ unknown elements. The log-likelihood function of the ECB regression model is expressed as

$$\begin{array}{ccc}\hfill \ell \left(\mathbf{\vartheta }\right)& =& nlog\left(log(\tau )\right)+nlog\left(log\left(\beta /(1-\beta )\right)\right)+log\left(\beta /(1-\beta )\right){\displaystyle \sum _{i=1}^n}{y}_i\hfill \\ & & -{\displaystyle \sum _{i=1}^n}log\left(log\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)}{\left(2\beta -1\right)}\right)\right)-{\displaystyle \sum _{i=1}^n}log\left({\left(\beta /(1-\beta )\right)}^{y_i}-1\right)\hfill \\ & & +log(\tau ){\displaystyle \sum _{i=1}^n}log\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{y_i}-1\right)}{\left(2\beta -1\right)}\right){\left(log\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)}{\left(2\beta -1\right)}\right)\right)}^{-1}.\hfill \end{array}$$

(17)

Thus, the first derivatives according to the model coefficients established in (17) are given by

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}\ell (\mathbf{\vartheta })}{\mathrm{d}{\beta }_r}& =& -{\displaystyle \sum _{i=1}^n}\frac{{\left(\beta /(1-\beta )\right)}^{{\mu }_i}{\mu }_i\left(1-{\mu }_i\right)x_{ir}log\left(\beta /(1-\beta )\right)}{\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)log\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)}{\left(2\beta -1\right)}\right)},r\in \{0,1,\cdots ,p\},\hfill \\ & & -{\displaystyle \sum _{i=1}^n}\frac{{\left(\beta /(1-\beta )\right)}^{{\mu }_i}{\mu }_i\left(1-{\mu }_i\right)x_{ir}log(\tau )log\left(\beta /(1-\beta )\right)}{{log}^2\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)}{\left(2\beta -1\right)}\right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)}log\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{y_i}-1\right)}{\left(2\beta -1\right)}\right),\hfill \\ \hfill \frac{\mathrm{d}\ell (\mathbf{\vartheta })}{\mathrm{d}\beta }& =& \frac{n}{\beta \left(1-\beta \right)log\left(\beta /(1-\beta )\right)}+\frac{1}{\beta \left(1-\beta \right)}{\displaystyle \sum _{i=1}^n}{y}_i-{\displaystyle \sum _{i=1}^n}\frac{y_i{\left(\beta /(1-\beta )\right)}^{y_i}}{\beta \left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{y_i}-1\right)}\hfill \\ & & -{\displaystyle \sum _{i=1}^n}\frac{{\left(\beta /(1-\beta )\right)}^{{\mu }_i}\left(\left(2\beta -1\right){\mu }_i-\beta \right)+\beta }{\beta \left(1-\beta \right)\left(2\beta -1\right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)log\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)}{\left(2\beta -1\right)}\right)}\hfill \\ & & +2log(\tau ){\displaystyle \sum _{i=1}^n}\frac{\left(\left(\beta \left(y_i-1/2\right)-\frac{y_i}{2}\right){\left(\beta /(1-\beta )\right)}^{y_i}+\frac{\beta }{2}\right)}{\beta \left(2\beta -1\right)\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{y_i}-1\right)log\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)}{\left(2\beta -1\right)}\right)}\hfill \\ & & -2log(\tau ){\displaystyle \sum _{i=1}^n}\frac{\left(\left(\beta \left({\mu }_i-1/2\right)-\frac{{\mu }_i}{2}\right){\left(\beta /(1-\beta )\right)}^{{\mu }_i}+\frac{\beta }{2}\right)log\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{y_i}-1\right)}{\left(2\beta -1\right)}\right)}{\beta \left(2\beta -1\right)\left(1-\beta \right){log}^2\left(\frac{\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)}{\left(2\beta -1\right)}\right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)}.\hfill \end{array}$$

Once the derivatives are set to zero, the resulting equations should be solved using numerical methods since they are nonlinear functions. The formula presented in (17) can be directly maximized utilizing the R software to obtain these solutions. It is important to note that the asymptotic distribution of $(\widehat{\mathbf{\vartheta }}-\mathbf{\vartheta })$ is a multivariate normal ${\mathrm{N}}_{p+2}(\mathbf{0},{\mathcal{I}}^{-1}(\mathbf{\vartheta }))$, where $\mathcal{I}(\mathbf{\vartheta })$ is the Fisher information matrix. However, for numerical applications, the $(p+2)\times (p+2)$ observed information matrix is often used instead of $\mathcal{I}(\mathbf{\vartheta })$ to make inference about the corresponding parameters.

4.4. Residual Analysis for Model Fitting

Validating a fitted regression model involves analyzing residuals. The randomized quantile (RQ) residual is a preferred method for this validation. The RQ residual, proposed in [32], is defined as

$${\widehat{r}}_i={\mathsf{\Phi }}^{-1}(G(y_i,\widehat{\beta },{\widehat{\mu }}_i,\tau )),i\in \{1,\cdots ,n\},$$

where $G(y,\beta ,\mu ,\tau )$ is the CDF of the CEB model and ${\widehat{\mu }}_i$ is defined in (16). If the fitted regression is valid, the RQ quantile should be standard normal distributed.

5. Simulation Studies

In this section, we discuss simulation studies for point estimation of the CB parameter and estimation of the ECB QR parameters.

5.1. Computational Framework

The experiments were conducted on a computer with Windows 10 for 64 bits, 8 Gigabytes of RAM, and an Intel Core i7-4510U 2-2.60 GHz processor. The optimization of the log-likelihood function can be performed directly by software such as R or Matlab.

In the case of using R, two libraries, named ConstrOptim and maxLik, can be employed to conduct the experimental results of the simulation studies and real data analysis when obtaining the point estimates of the corresponding parameters. The functions of these libraries numerically provide the parameter estimates and the associated estimated standard errors based on diverse algorithms. The Nelder–Mead optimization routine can be utilized to obtain these estimates. Such libraries may only be employed in specific contexts, which is their limitation. Our computations were carried out by the maxLik and goftest functions of the R software.

For more details about R packages related to QR, the interested reader is referred to a package named unitquantreg, which was developed in the setting stats::lm package [33]. The unitquantreg() function of this package is flexible and allows us, by means of the ns() function, to utilize regression splines. This function is available in the splines package. Moreover, it is possible to employ the gam() function of an R package named mgcv.

The unitquantreg package can be secured from https://github.com/AndrMenezes/unitquantreg (accessed on 5 April 2023) and installed using

• devtools::install_github(“AndrMenezes/unitquantreg”).

The unitquantreg package contains 15 distributions, and our new distribution can be added to this package.

5.2. Point Estimation

For the simulations of point estimation, we change n in $\{25,50,100,200,500,1000\}$ with $M=1000$ Monte Carlo replicates from the CB distribution, with $\beta \in \{0.25,0.75\}$. The empirical mean, bias, and mean square error (MSE) of the estimators, given by $\mathrm{Mean}(\widehat{\beta })=(1/M){\sum }_{i=1}^M{\widehat{\beta }}_i$, $\mathrm{Bias}(\widehat{\beta })=(1/M){\sum }_{i=1}^M(\beta -{\widehat{\beta }}_i)$, and $\mathrm{MSE}(\widehat{\beta })=(1/M){\sum }_{i=1}^M{(\beta -{\widehat{\beta }}_i)}^2$, were used to calculate statistical indicators. The results of the estimators are close to the true values of the parameters, and consistent, with MSE and bias decreasing as n increases. The moment estimate was used as the initial value for the ML estimates. The results are shown in Figure 5.

Simulation studies also were conducted to evaluate the behavior of 95% confidence intervals using the ML method. The coverage length (CL) and coverage probability (CP) were empirically calculated to assess this behavior. The CL and CP are defined by

$$\begin{array}{ccc}\hfill \mathrm{CL}(\widehat{\beta })& =& \frac{1}{M}\sum _{i=1}^{M}3.92s({\widehat{\beta }}_i),\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill \mathrm{CP}(\widehat{\beta })& =& \frac{1}{M}\sum _{i=1}^M𝟙({\widehat{\beta }}_i-1.96s({\widehat{\beta }}_i)<\beta <{\widehat{\beta }}_i+1.96s({\widehat{\beta }}_i)),\hfill \end{array}$$

(19)

respectively. In the expression given in (18), $s(\widehat{{\beta }_i})$ is the estimated standard error of the ML estimator. In the formula stated in (19), $𝟙$ is the indicator function. The simulation results for the empirical CLs and CPs are presented in Figure 6. It is evident from the figure that the empirical CLs decrease as the sample size increases, while the empirical CPs remain close to 0.95.

5.3. Quantile Regression Model

Next, we conduct simulations to evaluate the ML estimators’ performance for the ECB QR parameters. We measure their performance by examining the bias, MSE, CL, and CP across different sample sizes (n), known values of $\tau$, true values of $\beta$, and various covariate values. To generate the values of the unit response, we use a structure stated as

$$y_i=log\left(1+\frac{\left(2\beta -1\right)u_i^{1/{\alpha }_i}}{1-\beta }\right){\left(log\left(\beta /(1-\beta )\right)\right)}^{-1},i\in \{1,\cdots ,n\},$$

where $u_i\sim \mathrm{U}(0,1)$ and ${\alpha }_i=log(\tau {\left(log\left(\left(1-\beta \right)\left({\left(\beta /(1-\beta )\right)}^{{\mu }_i}-1\right)/\left(2\beta -1\right)\right)\right)}^{-1})$.

For this simulation study, we use $M=1000$ Monte Carlo replicates and varying sample sizes from $n=25$ to $n=1000$, with step sizes of $25,50,100,250,500$, and 1000. We also vary the values of $\tau$ and $\beta$, with $\tau$ taking on the values of $0.25$, $0.50$, and $0.75$ and $\beta$ taking on the values of $0.25$, $0.5$, and $0.75$. The regression structure used is formulated as

$$\mathrm{logit}({\mu }_i)={\beta }_0+{\beta }_1{z}_{i1},i\in \{1,\cdots ,n\},$$

with ${\beta }_0=-0.5$, ${\beta }_1=2$ and $z_{i1}\sim \mathrm{Bernoulli}(0.5)$.

The simulation results for the ECB QR model are presented in

Table 3. Empirical bias, MSE, $95\%$ CP, and $95\%$ CL for the ECB QR model with $\beta =0.25$ and $\tau =0.25$ for the indicated n and parameter with simulated data.

n	Bias(${\widehat{\mathit{\beta }}}_0$)	Bias(${\widehat{\mathit{\beta }}}_1$)	Bias($\widehat{\mathit{\beta }}$)	MSE(${\widehat{\mathit{\beta }}}_0$)	MSE(${\widehat{\mathit{\beta }}}_1$)	MSE($\widehat{\mathit{\beta }}$)	CL(${\widehat{\mathit{\beta }}}_0$)	CL(${\widehat{\mathit{\beta }}}_1$)	CL($\widehat{\mathit{\beta }}$)	CP(${\widehat{\mathit{\beta }}}_0$)	CP(${\widehat{\mathit{\beta }}}_1$)	CP($\widehat{\mathit{\beta }}$)
25	$0.0506$	$-0.0241$	$-0.0461$	$0.1578$	$0.2368$	$0.0524$	$1.3731$	$1.7408$	$0.8412$	$0.8880$	$0.9110$	$0.6430$
50	$0.0333$	$-0.0053$	$0.0057$	$0.0804$	$0.1169$	$0.0417$	$1.0516$	$1.3099$	$0.7733$	$0.9290$	$0.9490$	$0.7788$
100	$0.0120$	$-0.0009$	$-0.0064$	$0.0371$	$0.0593$	$0.0219$	$0.7197$	$0.9227$	$0.5588$	$0.9250$	$0.9430$	$0.8340$
200	$0.0078$	$-0.0008$	$-0.0013$	$0.0196$	$0.0298$	$0.0127$	$0.5310$	$0.6634$	$0.4251$	$0.9390$	$0.9360$	$0.8830$
500	$0.0004$	$0.0029$	$0.0001$	$0.0076$	$0.0122$	$0.0051$	$0.3351$	$0.4210$	$0.2758$	$0.9460$	$0.9510$	$0.9079$
1000	$-0.0010$	$-0.0001$	$0.0003$	$0.0035$	$0.0052$	$0.0027$	$0.2318$	$0.2967$	$0.1936$	$0.9500$	$0.9630$	$0.9240$

,

Table 4. Empirical bias, MSE, $95\%$ CP, and $95\%$ CL for the ECB QR model with $\beta =0.75$ and $\tau =0.25$ for the indicated n and parameter with simulated data.

n	Bias(${\widehat{\mathit{\beta }}}_0$)	Bias(${\widehat{\mathit{\beta }}}_1$)	Bias($\widehat{\mathit{\beta }}$)	MSE(${\widehat{\mathit{\beta }}}_0$)	MSE(${\widehat{\mathit{\beta }}}_1$)	MSE($\widehat{\mathit{\beta }}$)	CL(${\widehat{\mathit{\beta }}}_0$)	CL(${\widehat{\mathit{\beta }}}_1$)	CL($\widehat{\mathit{\beta }}$)	CP(${\widehat{\mathit{\beta }}}_0$)	CP(${\widehat{\mathit{\beta }}}_1$)	CP($\widehat{\mathit{\beta }}$)
25	$0.0227$	$0.0017$	$-0.2301$	$0.2340$	$0.3634$	$0.1553$	$1.7324$	$2.1644$	$1.1733$	$0.9010$	$0.9190$	$0.6650$
50	$0.0296$	$-0.0311$	$-0.1406$	$0.1236$	$0.1746$	$0.0915$	$1.3305$	$1.5991$	$0.9841$	$0.9180$	$0.9360$	$0.7730$
100	$0.0171$	$-0.0141$	$-0.0833$	$0.0748$	$0.0979$	$0.0531$	$0.9991$	$1.1703$	$0.7680$	$0.9180$	$0.9320$	$0.8138$
200	$0.0027$	$-0.0024$	$-0.0442$	$0.0335$	$0.0463$	$0.0249$	$0.6950$	$0.8258$	$0.5717$	$0.9500$	$0.9440$	$0.8940$
500	$0.0046$	$-0.0025$	$-0.0206$	$0.0131$	$0.0184$	$0.0100$	$0.4437$	$0.5251$	$0.3749$	$0.9430$	$0.9500$	$0.9259$
1000	$0.0014$	$-0.0020$	$-0.0120$	$0.0069$	$0.0097$	$0.0050$	$0.3280$	$0.3790$	$0.2774$	$0.9520$	$0.9420$	$0.9390$

,

Table 5. Empirical bias, MSE, $95\%$ CP, and $95\%$ CL for the ECB QR model with $\beta =0.50$ and $\tau =0.25$ for the indicated n and parameter with simulated data.

n	Bias(${\widehat{\mathit{\beta }}}_0$)	Bias(${\widehat{\mathit{\beta }}}_1$)	MSE(${\widehat{\mathit{\beta }}}_0$)	MSE(${\widehat{\mathit{\beta }}}_1$)	CL(${\widehat{\mathit{\beta }}}_0$)	CL(${\widehat{\mathit{\beta }}}_1$)	CP(${\widehat{\mathit{\beta }}}_0$)	CP(${\widehat{\mathit{\beta }}}_1$)
25	$0.0496$	$0.0033$	$0.0901$	$0.1766$	$1.1570$	$1.6256$	$0.9520$	$0.9480$
50	$0.0196$	$0.0110$	$0.0432$	$0.0875$	$0.7888$	$1.1547$	$0.9450$	$0.9500$
100	$0.0055$	$0.0012$	$0.0226$	$0.0465$	$0.5630$	$0.8153$	$0.9320$	$0.9410$
200	$0.0063$	$0.0032$	$0.0117$	$0.0214$	$0.4475$	$0.5820$	$0.9620$	$0.9540$
500	$0.0012$	$0.0001$	$0.0047$	$0.0088$	$0.2651$	$0.3638$	$0.9580$	$0.9470$
1000	$0.0022$	$0.0007$	$0.0025$	$0.0042$	$0.1919$	$0.2579$	$0.9420$	$0.9470$

,

Table 6. Empirical bias, MSE, $95\%$ CP, and $95\%$ CL for the ECB QR model with $\beta =0.50$ and $\tau =0.75$ for the indicated n and parameter with simulated data.

n	Bias(${\widehat{\mathit{\beta }}}_0$)	Bias(${\widehat{\mathit{\beta }}}_1$)	MSE(${\widehat{\mathit{\beta }}}_0$)	MSE(${\widehat{\mathit{\beta }}}_1$)	CL(${\widehat{\mathit{\beta }}}_0$)	CL(${\widehat{\mathit{\beta }}}_1$)	CP(${\widehat{\mathit{\beta }}}_0$)	CP(${\widehat{\mathit{\beta }}}_1$)
25	$0.0466$	$-0.0147$	$0.1133$	$0.1821$	$1.2587$	$1.6439$	$0.9250$	$0.9450$
50	$0.0220$	$-0.0064$	$0.0443$	$0.0876$	$0.8032$	$1.1513$	$0.9470$	$0.9470$
100	$0.0146$	$0.0064$	$0.0204$	$0.0418$	$0.5626$	$0.8150$	$0.9520$	$0.9580$
200	$0.0085$	$-0.0130$	$0.0120$	$0.0212$	$0.4259$	$0.5761$	$0.9480$	$0.9510$
500	$0.0053$	$-0.0041$	$0.0045$	$0.0087$	$0.2624$	$0.3636$	$0.9460$	$0.9480$
1000	$-0.0002$	$0.0004$	$0.0024$	$0.0045$	$0.1901$	$0.2576$	$0.9430$	$0.9420$

,

Table 7. Empirical bias, MSE, $95\%$ CP, and $95\%$ CL for the ECB QR model with $\beta =0.25$ and $\tau =0.50$ for the indicated n and parameter with simulated data.

n	Bias(${\widehat{\mathit{\beta }}}_0$)	Bias(${\widehat{\mathit{\beta }}}_1$)	Bias($\widehat{\mathit{\beta }}$)	MSE(${\widehat{\mathit{\beta }}}_0$)	MSE(${\widehat{\mathit{\beta }}}_1$)	MSE($\widehat{\mathit{\beta }}$)	CL(${\widehat{\mathit{\beta }}}_0$)	CL(${\widehat{\mathit{\beta }}}_1$)	CL($\widehat{\mathit{\beta }}$)	CP(${\widehat{\mathit{\beta }}}_0$)	CP(${\widehat{\mathit{\beta }}}_1$)	CP($\widehat{\mathit{\beta }}$)
25	$-0.0429$	$0.0319$	$-0.0166$	$0.2150$	$0.2637$	$0.0531$	$1.7293$	$1.8910$	$0.8087$	$0.9210$	$0.9260$	$0.6940$
50	$-0.0056$	$-0.0054$	$-0.0086$	$0.0798$	$0.1155$	$0.0305$	$1.1040$	$1.3311$	$0.6289$	$0.9380$	$0.9500$	$0.8060$
100	$-0.0070$	$0.0160$	$-0.0053$	$0.0431$	$0.0604$	$0.0147$	$0.8228$	$0.9506$	$0.4911$	$0.9550$	$0.9460$	$0.8700$
200	$-0.0007$	-$0.0009$	$-0.0023$	$0.0229$	$0.0310$	$0.0091$	$0.6091$	$0.6801$	$0.3655$	$0.9520$	$0.9490$	$0.8950$
500	$-0.0017$	-$0.0041$	$-0.0003$	$0.0096$	$0.0118$	$0.0038$	$0.3913$	$0.4336$	$0.2382$	$0.9560$	$0.9540$	$0.9300$
1000	$0.0008$	-$0.0013$	$0.0020$	$0.0049$	$0.0059$	$0.0019$	$0.2785$	$0.3076$	$0.1710$	$0.9590$	$0.9450$	$0.9420$

and

Table 8. Empirical bias, MSE, $95\%$ CP, and $95\%$ CL for the ECB QR model with $\beta =0.75$ and $\tau =0.50$ for the indicated n and parameter with simulated data.

n	Bias(${\widehat{\mathit{\beta }}}_0$)	Bias(${\widehat{\mathit{\beta }}}_1$)	Bias($\widehat{\mathit{\beta }}$)	MSE(${\widehat{\mathit{\beta }}}_0$)	MSE(${\widehat{\mathit{\beta }}}_1$)	MSE($\widehat{\mathit{\beta }}$)	CL(${\widehat{\mathit{\beta }}}_0$)	CL(${\widehat{\mathit{\beta }}}_1$)	CL($\widehat{\mathit{\beta }}$)	CP(${\widehat{\mathit{\beta }}}_0$)	CP(${\widehat{\mathit{\beta }}}_1$)	CP($\widehat{\mathit{\beta }}$)
25	$-0.0501$	$0.0264$	$-0.1891$	$0.2898$	$0.3932$	$0.1172$	$1.9476$	$2.3130$	$1.0495$	$0.9120$	$0.9290$	$0.7404$
50	$-0.0283$	$0.0257$	$-0.1072$	$0.1507$	$0.2058$	$0.0638$	$1.4803$	$1.6950$	$0.8336$	$0.9260$	$0.9320$	$0.8029$
100	$-0.0165$	$0.0129$	$-0.0566$	$0.0955$	$0.1078$	$0.0322$	$1.1702$	$1.2454$	$0.6562$	$0.9300$	$0.9400$	$0.8805$
200	$-0.0105$	$0.0059$	$-0.0321$	$0.0540$	$0.0589$	$0.0169$	$0.8995$	$0.9262$	$0.4952$	$0.9440$	$0.9360$	$0.8958$
500	$-0.0007$	$-0.0005$	$-0.0150$	$0.0188$	$0.0209$	$0.0063$	$0.5222$	$0.5580$	$0.3107$	$0.9370$	$0.9470$	$0.9333$
1000	$-0.0013$	$0.0012$	$-0.0056$	$0.0090$	$0.0101$	$0.0028$	$0.3792$	$0.4032$	$0.2310$	$0.9370$	$0.9500$	$0.9306$

. It is observed that the empirical CLs exhibit a decreasing trend as the sample size increases, while the empirical CPs approach 0.95. The biases are negligible, and all MSEs approach zero.

6. Applications in Science and Education

This section provides two real data illustrations to show potential applications of the CB distribution and ECB QR model. All datasets were collected via: OECD.Stat with link https://stats.oecd.org/ (accessed on 1 May 2023).

The OECD.Stat provides data on various development themes, including demography, education and training, health, labor, social protection and well-being, finance, as well as agriculture and fisheries, for both OECD and nonmember economies. Among these, we focus on two datasets. The first dataset includes the long-term unemployment rate of various countries, while the second dataset comprises the educational attainment proportion and better life index (BLI) measurements of different countries. We analyze the relationship between the educational attainment proportion and BLI using our QR model. All the datasets employed in this study are publicly available and can be secured directly from: https://stats.oecd.org/index.aspx?DataSetCode=BLI2017 (accessed on 1 May 2023).

6.1. Application I

We analyze the proportion of long-term unemployment in different countries, defined as the ratio of the number of individuals who have been unemployed for one year or more to the total labor force (employed and unemployed individuals). We compare the performance of four one-parameter distributions, namely CB, log-Bilal [12], unit Lindley [14], and Topp–Leone (TL) [34], using indicators such as the estimated log-likelihood ($\widehat{\ell }$), AIC (Akaike information criterion), and BIC (Bayesian information criterion), as well as Kolmogorov–Smirnov (KS), Cramer–von Mises (CM), and Anderson–Darling (AD) goodness-of-fit statistics.

We report the results of the data analysis utilizing ML estimation in

Table 9. ML estimate and standard errors, in parentheses, as well as $\widehat{\ell }$ and p-value of the goodness-of-fit test, in brackets, with unemployment data.

Distribution	$\widehat{\mathit{\beta }}$	$\widehat{\mathit{\ell }}$	AIC	BIC	AD	CM	KS
CB	$0.261\times {10}^{-13}$$(0.13\times {10}^{-12})$	$92.8248$	$-183.6497$	$-182.0122$	$0.6018$	$0.1058$	$0.1332\left[0.5097\right]$
Log-Bilal	$4.8700(0.5694)$	$74.2699$	$-146.5399$	$-144.9023$	$5.7487$	$1.1128$	$0.3077\left[0.0015\right]$
TL	$0.2996(0.0486)$	$68.0997$	$-134.1995$	$-132.5619$	$35.2670$	$7.0160$	$0.7110[<0.001]$
Unit-Lindley	$29.7750(4.8431)$	$92.2610$	$-182.5221$	$-180.8845$	$0.7448$	$0.1317$	$0.1532\left[0.3344\right]$

, which shows that the CB distribution is the best model according to the lowest values of AIC, BIC, AD, CM, and KS statistics. Plots of fitted PDF and CDF for the unemployment data are provided in Figure 7.

6.2. Application II

Now, we focus on the modeling ability of the ECB QR and compare its results with those of the KW [17] and LEEG [19] QRs. The CDFs of the KW and LEEG distributions are, respectively, given by

$$G_{\mathrm{KW}}(y;\beta ,\mu ,\tau )=1-{\left(1-y^{\beta }\right)}^{\frac{log(1-\tau )}{log(1-{\mu }^{\beta })}}y\in (0,1),$$

$$G_{\mathrm{LEEG}}(y;\beta ,\mu ,\tau )=\frac{\tau y^{\beta }\left({\mu }^{\beta }-1\right)}{\tau \left({\mu }^{\beta }-y^{\beta }\right)+{\mu }^{\beta }\left(y^{\beta }-1\right)},y\in (0,1),$$

where $\mu \in (0,1)$ is the quantile parameter, $\beta >0$ is a parameter related to the shape of the distribution, and $\tau$ is the known quantile level.

We investigate the relationship between the conditional quantile of educational attainment and the independent variables, including homicide rate (HR), dwellings without basic facilities (DWBF), and labor market insecurity (LMI) values of countries. The unit-dependent variable employed in this study is the proportion of educational attainment. This dataset was previously analyzed in [12] using the log-Bilal distribution for unit mean response regression modeling. The QR based on ${\mu }_i$ is given by

$${\mu }_i=\frac{{\mathrm{e}}^{{\beta }_0+{\beta }_1{\mathrm{HR}}_i+{\beta }_2{\mathrm{DWBF}}_i+{\beta }_3{\mathrm{LMI}}_i}}{1+{\mathrm{e}}^{{\beta }_0+{\beta }_1{\mathrm{HR}}_i+{\beta }_2{\mathrm{DWBF}}_i+{\beta }_3{\mathrm{LMI}}_i}},i\in \{1,\cdots ,38\},$$

where ${\mu }_i$ is the quantile for all models.

We obtain the QR results for different quantile levels such as $\tau =0.25$, $\tau =0.5$, and $\tau =0.75$, separately.

Table 10. Results of fitted regressions with quantile $\tau =0.25$ and model selection criteria for educational data.

Parameter	ECB			KW			LEEG
	Estimate	SE	p-Value	Estimate	SE	p-Value	Estimate	SE	p-Value
${\beta }_0$	1.9160	0.3172	<0.0001	1.3473	0.1730	<0.0001	1.4199	0.1599	<0.0001
${\beta }_1$	−0.0914	0.0241	<0.0001	−0.0511	0.0130	<0.0001	−0.05588	0.0118	<0.0001
${\beta }_2$	11.1800	2.3920	<0.0001	2.0563	1.1999	0.0866	2.3612	1.3068	0.078
${\beta }_3$	−21.2000	6.1390	<0.001	−7.1977	1.1596	<0.0001	−8.1841	1.7661	<0.0001
$\beta$	0.00058	0.0009	0.5367	7.0769	1.2104	<0.0001	9.1738	1.9582	<0.0001
$\widehat{\ell }$	35.5794			33.5838			32.3569
AIC	−61.1587			−57.1676			−54.7139
BIC	−52.9708			−48.9797			−46.5259

,

Table 11. Results of fitted regressions with quantile $\tau =0.50$ and model selection criteria for educational data.

Parameter	ECB			KW			LEEG
	Estimate	SE	p-Value	Estimate	SE	p-Value	Estimate	SE	p-Value
${\beta }_0$	2.3830	1.7410	<0.0001	1.9595	0.1795	<0.0001	1.9907	0.1903	<0.0001
${\beta }_1$	−0.0942	0.0207	<0.0001	−0.0639	0.0135	<0.0001	−0.0671	0.0134	<0.0001
${\beta }_2$	10.9300	2.2890	<0.0001	2.3783	1.2035	0.0481	2.8145	1.4914	0.0591
${\beta }_3$	−21.0700	1.8860	<0.0001	−8.9696	1.3384	<0.0001	−9.8119	2.0626	<0.0001
$\beta$	0.00097	0.0015	0.513	7.2166	1.2191	<0.0001	9.1091	1.9682	<0.0001
$\widehat{\ell }$	34.6717			33.8642			32.2079
AIC	−59.3433			−57.7286			−54.4158
BIC	−51.1553			−49.5406			−46.2279

and

Table 12. Results of fitted regressions with quantile $\tau =0.75$ and model selection criteria for educational data.

Parameters	ECB			KW			LEEG
	Estimate	SE	p-Value	Estimate	SE	p-Value	Estimate	SE	p-Value
${\beta }_0$	3.0897	0.2262	<0.0001	2.5204	0.2044	<0.0001	2.7506	0.2707	<0.0001
${\beta }_1$	−0.0978	0.0229	<0.0001	−0.0772	0.0133	<0.0001	−0.0834	0.0157	<0.0001
${\beta }_2$	9.7876	1.8659	<0.0001	2.7295	1.5489	0.0708	3.5070	1.8068	0.0523
${\beta }_3$	−19.9055	1.4959	<0.0001	−10.0828	1.6659	<0.0001	−12.1760	2.6050	<0.0001
$\beta$	0.0022	0.0029	0.4530	8.2613	1.3396	<0.0001	9.0614	2.0043	<0.0001
$\widehat{\ell }$	33.8923			33.6323			32.0757
AIC	−57.7847			−57.2647			−54.1515
BIC	−49.5967			−49.0768			−45.9635

supply the results of the QR analysis for the different quantile levels mentioned. The findings from these tables indicate that when inferring about the ECB QR, all independent variables are statistically significant at any usual significance level based on all quantile levels.

The signs of ${\beta }_0$ and ${\beta }_2$ are positive, whereas the signs of ${\beta }_1$ and ${\beta }_3$ are negative. Hence, all independent variables have effects on the educational attainment of the countries. Since a direct positive correlation between the educational attainment and DWBF exists, as the educational attainment of countries increases, the DWBF increases as well. Moreover, there are negative relationships between educational attainment and the HR and LMI index, indicating that as the educational attainment of countries increases, the HR and LMI index decreases. Furthermore, the ECB QR model has minimal AIC/BIC with maximal log-likelihood values compared with other regression models at all quantile levels. Consequently, we conclude that the best model is the proposed regression structure.

Figure 8, Figure 9 and Figure 10 display the QQ (empirical quantile versus theoretical quantile) plot of the RQ residual for all QRs based on all quantile levels. These figures show that explaining the unit-dependent variable by independent variables is remarkable, and there is a good fit for the ECB QR model.

7. Conclusions

This article presented new properties of the continuous Bernoulli distribution and introduced its quantile regression model with an exponentiated structure. Our quantile regression model incorporated independent variables that were associated with the quantile of the dependent variable using a link function related to a logit structure. To assess the performance of the estimators, we employed the Monte Carlo method under various settings. Our results indicated that the estimators work well in statistical terms, and additionally, the coverage probabilities approach the nominal value as the size of the sample increases.

The adequacy of the continuous Bernoulli distribution, exponentiated continuous Bernoulli quantile regression model, and other regression structures was compared by analyzing two real datasets based on the better life index of the countries. The comparison was carried out based on statistical indicators. Results indicated that the continuous Bernoulli distribution and its associated quantile regression model demonstrated a good coherence with the two sets based on real data, providing an alternative to other models presented in the literature on the topic.

Our quantile regression model yielded similar results to the unit mean regression models presented in [12,13]. Other related works [35,36,37] proposed the unit log-log, transmuted unit Rayleigh, and unit quantile regression models, respectively. These models analyzed the educational attainment values of OECD countries with a unit-dependent variable and reported a negative relationship between all responses and the independent variable of homicide rate, which is consistent with the findings of our study.

The examination we conducted showed that the continuous Bernoulli distribution can be very effective in analyzing data that are limited to the unit interval, whether or not independent variables are present. The use of both quantile and mean regression models can be easily carried out through the implementation of R codes by the authors of this article.

There are several avenues for future research that can build upon the present investigation. One possible direction is to incorporate multivariate, functional, temporal, and spatial structures, errors-in-variables, and partial least squares into the quantile regression framework, as well as studying their influence on diagnostics [38]. Other potential areas of investigation are the Tobit and Cobb–Douglas frameworks [39], which could be relevant to the topic of this study. Additionally, the analysis of censored observations and frailty models in the present context is another avenue for further research [40]. The authors are currently analyzing these issues and their findings are expected to be available in the future.