Laplace-Logistic Unit Distribution with Application in Dynamic and Regression Analysis

1. Introduction

Unit stochastic distributions, defined on the interval $(0,1)$, represent an important tool of contemporary probability theory. They are primarily used as stochastic models that can describe so-called proportional (percentage) variables and represent theoretical models that can successfully explain the behaviour of some real-world phenomena (see, as more contemporary, e.g. [1,2,3,4,5,6,7,8,9,10,11,12,13]). On the other hand, modeling with unit distributions differs from common stochastic modeling procedures, primarily due to the limitation of data within $(0,1)$ interval. Although the procedure for creating unit distributions can be given in a general form [14], the most common approach is based on continuous transformations of distributions defined on infinite intervals into a unit interval (as some recent results, see e.g. [15,16,17,18,19,20,21,22,23]). Nevertheless, most of these distributions are limited by a unique form of (a)symmetry and modality, "vanishing" data at the ends of the unit interval, etc. This is often inappropriate for modeling real phenomena with different characteristics, especially where data with pronounced "peaks" and "fat tails" appear.

To this end, proceeding from similar considerations as in Stojanović et al. [24], Stojanović et al. [25], a new unit distribution, named the Laplace-logistic unit (LLU) distribution, is described here. It is based on a general logistic mapping of the Laplace distribution to a unit interval, and as will be seen, has considerable flexibility and suitability for describing a variety of empirical distributions, from those with pronounced extremes to increasing, decreasing or bathtub-shaped ones. The definition and the key stochastic properties of the LLU distribution, regarding its limit properties, (a)symmetry and modality, are presented in the next Section 2. Thereafter, Section 3 considers the parameters estimation based on sample quantiles, the asymptotic properties of thus obtained estimators, as well as their numerical Monte Carlo study. An application of the LLU distribution in fitting some real-world data, primarily from the aspect of dynamic and regression analysis, is presented in Section 4, while Section 5 contains some concluding remarks.

2. The LLU Distribution

2.1. Definition and Main Properties

Let us first consider a random variable (RV) Y with a symmetric Laplace distribution, whose probability density function (PDF) is:

$$g(y;\theta )=\frac{\lambda }{2}exp\left(-\lambda \left|y\right|\right),$$

(1)

where $y\in \mathbb{R}$ and $\lambda >0$ is the scale parameter. For an arbitrary $\theta >0$, let us define the so-called general logistic map $x=\phi (y;\theta )={\left(1+exp\left(y\right)\right)}^{-1/\theta }$, which is the bijective, continuous transformation $\phi :(-\infty ,+\infty )\times (0,+\infty )\to (0,1)$, with limits:

$$\underset{y\to -\infty }{lim}\phi (y;\theta )=1^{-},\underset{y\to +\infty }{lim}\phi (y;\theta )=0^{+}.$$

Using Equation (1) and the inverse transformation $y={\phi }^{-1}(x;\theta )=ln(x^{-\theta }-1)$, the new RV $X=\phi (Y;\theta )$, defined on the unit interval $(0,1)$ is obtained. After some computation, for the PDF of the RV X one obtains:

$$\begin{array}{cc}\hfill f(x;\lambda ,\theta )& =g\left({\phi }^{-1}(x;\theta );\lambda \right)·\left|\frac{\partial {\phi }^{-1}(x;\theta )}{\partial x}\right|=\frac{\theta \lambda }{2x(1-x^{\theta })}exp\left(-\lambda |ln(x^{-\theta }-1)|\right)\hfill \\ \hfill & =\left\{\begin{array}{cc}{\displaystyle \frac{\theta \lambda x^{\theta \lambda -1}}{2{\left(1-x^{\theta }\right)}^{\lambda +1}},}\hfill & x\le 2^{-1/\theta }\hfill \\ {\displaystyle \frac{\theta \lambda {\left(1-x^{\theta }\right)}^{\lambda -1}}{2x^{\theta \lambda +1}},}\hfill & x>2^{-1/\theta }\hfill \end{array}\right.,\hfill \end{array}$$

(2)

where $x\in (0,1)$ and ${\int }_0^{1}f(x;\lambda ,\theta )\mathrm{d}x=1$. Thus, we say that the RV X, whose PDF is given by Equation (2), has a Laplace-logistic unit (LLU) distribution, with the parameters $\lambda ,\theta >0$, or abbr. $X:\mathcal{L}(\lambda ,\theta )$. Obviously, LLU distribution is a two-parameter distribution, where, in addition to the scale parameter $\lambda >0$, there is also a shape parameter $\theta >0$. Also, let us emphasize that symmetric Laplace density is expressed in terms of the absolute difference from the zero and it has pronounced “peaks” and tails more "fat" than, for instance, the Gaussian distribution. For these reasons, the RV $X:\mathcal{L}(\lambda ,\theta )$ will have similar properties, but it also has some other specificities. To describe them more completely, we first introduce some terms related to the limit behaviour of the density at the ends of the unit interval:

Definition 1. 

Let $X:\mathcal{L}(\lambda ,\theta )$ be the RV with the LLU distribution, whose PDF $f(x;\lambda ,\theta )$ is given by Equation (2). The PDF $f(x;\lambda ,\theta )$ is left-(right-)vanishing if is valid:

$$\underset{x\to 0^{+}}{lim}f(x;\lambda ,\theta )=0^{+}\left(\underset{x\to 1^{-}}{lim}f(x;\lambda ,\theta )=0^{+}\right).$$

Otherwise, the PDF $f(x;\lambda ,\theta )$ is left-(right-)tailed.

Now, some properties of the LLU distribution regarding its possible shapes and boundary characteristics can be given by the following proposition:

Theorem 1. 

The PDF $f(x;\lambda ,\theta )$ of the RV $X:\mathcal{L}(\lambda ,\theta )$ is continuous function, non-differentiable at the point $x_0=2^{-1/\theta }$, with the following properties:

$\left(i\right)$

When $\theta (1+2\lambda )<1$ and $\lambda >1$, it is decreasing.

$\left(ii\right)$

When $\theta (1-2\lambda )>1$ and $\theta \lambda >1$, it is increasing.

$\left(iii\right)$

When $\theta (1+2\lambda )\ge 1$, $\theta \lambda \le 1$ and $\lambda >1$, it is left-tailed and right-vanishing.

$\left(iv\right)$

When $\theta (1-2\lambda )\le 1$, $\theta \lambda >1$ and $\lambda \le 1$, it is right-tailed and left-vanishing.

$\left(v\right)$

When $\theta \lambda \le 1$ and $\lambda \le 1$, it is (both sides) tailed.

$\left(vi\right)$

When $\theta \lambda >1$ and $\lambda >1$, it is (both sides) vanishing.

Proof. 

For simplicity, let us denote the left and right branches of the PDF $f(x;\lambda ,\theta )$, respectively, as follows:

$${\displaystyle f_1(x;\lambda ,\theta ):=\frac{\theta \lambda x^{\theta \lambda -1}}{2{\left(1-x^{\theta }\right)}^{\lambda +1}},f_2(x;\lambda ,\theta ):=\frac{\theta \lambda {\left(1-x^{\theta }\right)}^{\lambda -1}}{2x^{\theta \lambda +1}}.}$$

Accordingly, it is obtained:

$$\underset{x\to 2^{-1/\theta }}{\lim }{f}_1(x;\lambda ,\theta )=\underset{x\to 2^{-1/\theta }}{\lim }{f}_2(x;\lambda ,\theta )=\theta \lambda 2^{1/\theta },$$

so the PDF $f(x;\lambda ,\theta )$ is indeed continuous at $x_0=2^{-1/\theta }$. Furthermore, after some computation, one obtains:

$$\underset{x\to 0^{+}}{\lim }{f}_1(x;\lambda ,\theta )=\left\{\begin{array}{cc}0,\hfill & \theta \lambda >1\hfill \\ 1/2,\hfill & \theta \lambda =1\hfill \\ +\infty ,\hfill & \theta \lambda <1\hfill \end{array}\right.,\underset{x\to 1^{-}}{\lim }{f}_2(x;\lambda ,\theta )=\left\{\begin{array}{cc}0,\hfill & \lambda >1\hfill \\ \theta /2,\hfill & \lambda =1\hfill \\ +\infty ,\hfill & \lambda <1\hfill \end{array}\right.,$$

(3)

as well as:

$$\begin{array}{cc}\hfill \frac{\partial f_1(x;\lambda ,\theta )}{\partial x}& =\frac{\theta \lambda x^{\theta \lambda -2}}{2{\left(1-x^{\theta }\right)}^{\lambda +2}}\left((\theta +1)x^{\theta }+\theta \lambda -1\right),\hfill \\ \hfill \frac{\partial f_2(x;\lambda ,\theta )}{\partial x}& =\frac{\theta \lambda {\left(1-x^{\theta }\right)}^{\lambda -2}}{2x^{\theta \lambda +2}}\left((\theta +1)x^{\theta }-\theta \lambda -1\right).\hfill \end{array}$$

(4)

Thus, the sign of the partial derivatives in Equation (4) depends on the sign of the functions:

$${\psi }_1(x;\lambda ,\theta )=(\theta +1)x^{\theta }+\theta \lambda -1,{\psi }_2(x;\lambda ,\theta )=(\theta +1)x^{\theta }-\theta \lambda -1.$$

Obviously, both of these functions are monotonically increasing on $x\in (0,1)$, and their values at the critical points $x\in \{0,2^{-1/\theta },1\}$ are, respectively,

$$\begin{array}{cc}\hfill & {\psi }_1(0;\lambda ,\theta )=\theta \lambda -1,{\psi }_1(2^{-1/\theta };\lambda ,\theta )=\frac{1}{2}\left(\theta (2\lambda +1)-1\right),\hfill \\ \hfill & {\psi }_2(2^{-1/\theta };\lambda ,\theta )=\frac{1}{2}\left(\theta (1-2\lambda )-1\right),{\psi }_2(1;\lambda ,\theta )=\theta (1-\lambda ).\hfill \end{array}$$

(5)

Then, by using Equations (3) and (5), all cases in the statement of the theorem are simply obtained. □

Remark 1. 

The different conditions for the shape and boundary behavior of the LLU distribution can be seen in Figure 1 (a), where the six different areas mentioned in the previous Theorem are shown. On the other hand, Figure 1 (b) shows typical, different shapes of the PDFs of this distribution. As can be easily seen, the LLU distribution takes very different forms, where in addition to the typical one with a "peak", which is similar to the Laplace distribution, it can have a decreasing, increasing or bathtub-shaped PDF. In that way, it has a significant flexibility, which is of particular importance in applications. Moreover, when $\lambda \approx 0$ one obtains:

$${\psi }_1(x;\lambda ,\theta ){|}_{x=2^{-1/\theta }}\approx {\psi }_2(x;\lambda ,\theta ){|}_{x=2^{-1/\theta }}\approx 0,$$

so in this asymptotic case the PDF $f(x;\lambda ,\theta )$ will be an approximately differentiable, with extreme value (i.e. mode) at the point $x_0=2^{-1/\theta }$. In any case, the modality of the LLU distribution will be explored in more detail in the next part of this section.

The PDF of the LLU distribution can also be used to obtain its moments and other moment-based features, as given below.

Theorem 2. 

The $r^{th}$ moment of the LLU of the distributed RV X can be expressed as:

$$\begin{array}{cc}\hfill {\mu }_r(\lambda ,\theta ):=E\left(X^r\right)& =\frac{\lambda }{2}\left(B_{1/2}\left(\frac{r}{\theta }+\lambda ,-\lambda \right)+B_{1/2}\left(\lambda ,\frac{r}{\theta }-\lambda \right)\right),\hfill \end{array}$$

(6)

where r is an integer and $B_z(a,b)={\int }_0^z{t}^{a-1}{(1-t)}^{b-1}\mathrm{d}t$ is an incomplete beta function. In addition, the following convergence holds:

$$\underset{\theta \to \infty }{lim}{\mu }_r(\lambda ,\theta )=1.$$

(7)

Proof. 

By definition of the LLU distribution, its $r^{th}$ moment can be computed as follows:

$${\mu }_r(\lambda ,\theta )={\int }_0^1{x}^{r}f(x;\lambda ,\theta )\mathrm{d}x={\mu }_r^{\left(1\right)}(\lambda ,\theta )+{\mu }_r^{\left(2\right)}(\lambda ,\theta ),$$

(8)

wherein is:

$$\begin{array}{cc}\hfill {\mu }_r^{\left(1\right)}(\lambda ,\theta )& :={\int }_0^{x_0}{x}^r{f}_1(x;\lambda ,\theta )\mathrm{d}x=\frac{\lambda \theta }{2}{\int }_0^{x_0}\frac{x^{r+\lambda \theta -1}\mathrm{d}x}{{(1-x^{\theta })}^{\lambda +1}}=\frac{\lambda }{2}{\int }_0^{1/2}{t}^{r/\theta +\lambda -1}{(1-t)}^{-\lambda -1}\mathrm{d}t\hfill \\ \hfill & =\frac{\lambda }{2}{B}_{1/2}\left(\frac{r}{\theta }+\lambda ,-\lambda \right),\hfill \\ \hfill {\mu }_r^{\left(2\right)}(\lambda ,\theta )& :={\int }_{x_0}^1{x}^r{f}_2(x;\lambda ,\theta )\mathrm{d}x=\frac{\lambda \theta }{2}{\int }_{x_0}^1\frac{{(1-x^{\theta })}^{\lambda -1}\mathrm{d}x}{x^{\lambda \theta +1-r}}=\frac{\lambda }{2}{\int }_{1/2}^1{t}^{r/\theta -\lambda -1}{(1-t)}^{\lambda -1}\mathrm{d}t\hfill \\ \hfill & ={\int }_0^{1/2}{z}^{\lambda -1}{(1-z)}^{r/\theta -\lambda -1}\mathrm{d}z=\frac{\lambda }{2}{B}_{1/2}\left(\lambda ,\frac{r}{\theta }-\lambda \right)),\hfill \end{array}$$

as well as $x_0=2^{-1/\theta }$, $t=x^{\theta }$, $z=1-t$, and $B_z(a,b)={\int }_0^z{t}^{a-1}{(1-t)}^{b-1}\mathrm{d}t$. According to this and Equation (8), the first part of the theorem, that is Equation (6) is proven. Finally, Equation (7) follows from equalities:

$$B_{1/2}\left(\lambda ,-\lambda \right)={\int }_0^{1/2}\frac{t^{\lambda -1}\mathrm{d}t}{{(1-t)}^{\lambda +1}}={\int }_0^1{u}^{\lambda -1}\mathrm{d}u=\frac{1}{\lambda },$$

where $u=t/(1-t)$. □

Remark 2. 

By using Equation (6), for the mean and the variance of the RV $X:\mathcal{L}(\lambda ,\theta )$ ones obtain, respectively,

$$\begin{array}{cc}\hfill E\left(X\right)& ={\mu }_1(\lambda ,\theta )=\frac{\lambda }{2}\left[B_{\frac{1}{2}}\left(\lambda ,\frac{1}{\theta }-\lambda \right)+B_{\frac{1}{2}}\left(\lambda +\frac{1}{\theta },-\lambda \right)\right],\hfill \\ \hfill \mathrm{Var}\left(X\right)& ={\mu }_2(\lambda ,\theta )-{\left({\mu }_1(\lambda ,\theta )\right)}^2=\frac{\lambda }{2}[B_{\frac{1}{2}}\left(\lambda ,\frac{2}{\theta }-\lambda \right)+B_{\frac{1}{2}}\left(\lambda +\frac{2}{\theta },-\lambda \right)\hfill \\ \hfill & -\frac{\lambda }{2}{\left(B_{\frac{1}{2}}\left(\lambda ,\frac{1}{\theta }-\lambda \right)+B_{\frac{1}{2}}\left(\lambda +\frac{1}{\theta },-\lambda \right)\right)}^2].\hfill \end{array}$$

Additionally, according to Equation (7), it follows:

$$\underset{\theta \to \infty }{lim}E\left(X\right)=1,\underset{\theta \to \infty }{lim}\mathrm{Var}\left(X\right)=0,$$

that is, in this limit case, $X\stackrel{as}{=}1$ holds, where "as" means "almost surely". Therefore, the RV X is then reduced to a unit constant, obtained by transforming the Laplace distribution, whose PDF is given by Equation (1), using the trivial logistic map $\phi (y;\theta )\equiv 1$. This can also be seen in Figure 2 below, where 3D plots of mean value and variance are shown, depending on the parameters $\lambda ,\theta >0$.

Similarly, the skewness coefficient and the kurtosis of the RV $X:\mathcal{L}(\lambda ,\theta )$ are, respectively,

$$\begin{array}{cc}\hfill S(\lambda ,\theta )& :=\frac{E{(X-{\mu }_1(\lambda ,\theta ))}^3}{{\left[Var\left(X\right)\right]}^{3/2}}={\displaystyle \frac{{\mu }_3(\lambda ,\theta )-3{\mu }_1(\lambda ,\theta )\mathrm{Var}\left(X\right)-{\left({\mu }_1(\lambda ,\theta )\right)}^3}{{\left[Var\left(X\right)\right]}^{3/2}}},\hfill \\ \hfill K(\lambda ,\theta )& =\frac{E{\left(X-{\mu }_1(\lambda ,\theta )\right)}^4}{{\left[Var\left(X\right)\right]}^2}\hfill \\ \hfill & ={\displaystyle \frac{{\mu }_4(\lambda ,\theta )-4{\mu }_3(\lambda ,\theta ){\mu }_1(\lambda ,\theta )+6{\mu }_2(\lambda ,\theta ){\left({\mu }_1(\lambda ,\theta )\right)}^2-3{\left({\mu }_1(\lambda ,\theta )\right)}^4}{{\left[Var\left(X\right)\right]}^2}}.\hfill \end{array}$$

Nevertheless, because of the complexity though, a more detailed procedure for calculating these coefficients will be omitted.

2.2. Cumulative, Hazard and Quantile Function

Using Equation (2), the cumulative distribution function (CDF) of the LLU distribution can be obtained as follows:

$$F(x;\lambda ,\theta )=P\{X<x\}={\int }_0^{x}f(t;\lambda ,\theta )\mathrm{d}t=\left\{\begin{array}{cc}{\displaystyle \frac{x^{\theta \lambda }}{2{\left(1-x^{\theta }\right)}^{\lambda }}},\hfill & x\le 2^{-1/\theta }\hfill \\ 1-{\displaystyle \frac{{\left(1-x^{\theta }\right)}^{\lambda }}{2x^{\theta \lambda }}},\hfill & x>2^{-1/\theta }\hfill \end{array}\right.,$$

(9)

where $x\in (0,1)$. Note that the function $F(x;\lambda ,\theta )$ is differentiable on $x\in (0,1)$, and well defined outside the unit interval, since it is valid:

$$\underset{x\to 0^{+}}{lim}F(x;\lambda ,\theta )=0^{+},\underset{x\to 1^{-}}{lim}F(x;\lambda ,\theta )=1^{-}.$$

Furthermore, according to Equations (2) and (8), the hazard rate function (HRF) can be obtained as follows:

$$H(x;\lambda ,\theta ):=\frac{f(x;\lambda ,\theta )}{1-F(x;\lambda ,\theta )}=\left\{\begin{array}{cc}{\displaystyle \frac{\theta \lambda x^{\theta \lambda -1}}{\left(1-x^{\theta }\right)\left(2{\left(1-x^{\theta }\right)}^{\lambda }-x^{\theta \lambda }\right)}},\hfill & x\le 2^{-1/\theta }\hfill \\ {\displaystyle \frac{\theta \lambda }{x(1-x^{\theta })}},\hfill & x>2^{-1/\theta }\hfill \end{array}\right..$$

(10)

The basic properties of the HRF $H(x;\lambda ,\theta )$ can be expressed by the following statement:

Theorem 3. 

Let $X:\mathcal{L}(\lambda ,\theta )$ and $H(x;\lambda ,\theta )$ is the HRF of the RV X, defined by Equation (10). Then, $H(x;\lambda ,\theta )$ is a continuous function, non-differentiable at the point $x_0=2^{-1/\theta }$, with the following properties:

$\left(i\right)$

When $\theta \lambda \le 1$ and $\theta <1$, it is (both sides) tailed with a local maxima at $x_0=2^{-1/\theta }$.

$\left(ii\right)$

When $\theta \lambda \le 1$ and $\theta \ge 1$, it is (both sides) tailed without maxima, that is, bathtub-shaped.

$\left(iii\right)$

When $\theta \lambda >1$ and $\theta <1$, it is left-vanishing with a local maxima at $x_0=2^{-1/\theta }$.

$\left(iv\right)$

When $\theta \lambda >1$ and $\theta \ge 1$, it is increasing.

Proof. 

Similarly as in the proof of Theorem 1, let us denote:

$$H_1(x;\lambda ,\theta ):={\displaystyle \frac{\theta \lambda x^{\theta \lambda -1}}{\left(1-x^{\theta }\right)\left(2{\left(1-x^{\theta }\right)}^{\lambda }-x^{\theta \lambda }\right)}},H_2(x;\lambda ,\theta ):={\displaystyle \frac{\theta \lambda }{x(1-x^{\theta })}},$$

according to which it follows:

$$\underset{x\to 0^{+}}{\lim }{H}_1(x;\lambda ,\theta )=\left\{\begin{array}{cc}0,\hfill & \theta \lambda >1\hfill \\ 1/2,\hfill & \theta \lambda =1\hfill \\ +\infty ,\hfill & \theta \lambda <1\hfill \end{array}\right.,\underset{x\to 1^{-}}{\lim }{H}_2(x;\lambda ,\theta )=+\infty ,$$

(11)

as well as:

$$\begin{array}{cc}\hfill \frac{\partial H_1(x;\lambda ,\theta )}{\partial x}& =\frac{\theta \lambda x^{\theta \lambda -2}\left(x^{\theta \lambda }\left(1-(\theta +1)x^{\theta }\right)+2(\theta +1)x^{\theta }{\left(1-x^{\theta }\right)}^{\lambda }+2(\theta \lambda -1){\left(1-x^{\theta }\right)}^{\lambda }\right)}{{\left(x^{\theta }-1\right)}^2{\left(x^{\theta \lambda }-2{\left(1-x^{\theta }\right)}^{\lambda }\right)}^2},\hfill \\ \hfill \frac{\partial H_2(x;\lambda ,\theta )}{\partial x}& =\frac{\theta \lambda \left((\theta +1)x^{\theta }-1\right)}{x^2{\left(x^{\theta }-1\right)}^2}.\hfill \end{array}$$

(12)

Thus, Equations (11) and (12) give the statement of the theorem. □

Remark 3. 

Plots of the CDF and HRF of the LLU distribution are given in Figure 3 below. It is known that the hazard (failure) rate represents the failure frequency of the designed system or component. Thus, the HRF usually increases, which means that the probability that the designed system or component will fail increases. In contrary, decreasing failure rate (DFR) describes the phenomenon where this probability decreases in some interval. As can be concluded from the previous theorem, both of these situations can be obtained from the LLU distribution, for certain values of its parameters. Moreover, in cases where HRF has a local maximum at $x_0\in (0,1)$, the "peaked" value can then be interpreted as the critical point of the system. Therefore, these properties of HRF give diverse possibilities of its practical application, which will be discussed later.

In the last part of this section, the so-called quantile function (QF) of the LLU distribution is considered, obtained as inverse function of its CDF:

$$Q(p;\lambda ,\theta ):=F^{-1}(p;\lambda ,\theta )=\left\{\begin{array}{cc}{\left(2p\right)}^{\frac{1}{\theta \lambda }}{\left({\left(2p\right)}^{1/\lambda }+1\right)}^{-\frac{1}{\theta }},\hfill & 0<p\le 1/2\hfill \\ {\left({\left(2(1-p)\right)}^{1/\lambda }+1\right)}^{-\frac{1}{\theta }},\hfill & 1/2<p<1\hfill \end{array}\right.,$$

(13)

where $p=F(x;\lambda ,\theta )$. The QF is a useful tool for obtaining some other properties of the LLU distribution, primarily related to its modality and (a)symmetry.

Theorem 4. 

Let $X:\mathcal{L}(\lambda ,\theta )$ be the LLU distributed RV, whose QF is given by Equation (13). Then, the following statements hold:

$\left(i\right)$

The RV X is symmetrically distributed if and only if $\theta =1$. Otherwise, X is positively asymmetric when $0<\theta <1$, and negatively asymmetric when $\theta >1$.

$\left(ii\right)$

The RV X is unimodal, with the mode $x_0=2^{-1/\theta }$, if and only if $\theta (1+2\lambda )\ge 1$ and $\theta (1-2\lambda )\le 1$.

Proof. $\left(i\right)$ By substituting the quantile $p=1/2$ in the QF $Q(p;\lambda ,\theta )$, it is obtained the median $m=Q(1/2;\lambda ,\theta )=2^{-1/\theta }$. Thus, the RV X is symmetrically distributed if and only if the median is equal to 1/2, that is, when $\theta =1$. The positive and negative asymmetry conditions are also easily obtained, by solving the inequalities $Q(1/2;\lambda ,\theta )<1/2$ and $Q(1/2;\lambda ,\theta )>1/2$ , respectively.

$\left(ii\right)$ Using the rule for the derivative of the inverse function, for the derivatives up to the second order of the QF $Q(p;\lambda ,\theta )$ one obtains:

$$\begin{array}{cc}\hfill \frac{\partial Q(p;\lambda ,\theta )}{\partial p}& =\frac{1}{\partial F(x;\lambda ,\theta )/\partial x}=\frac{1}{f(x;\lambda ,\theta )},\hfill \\ \hfill \frac{\partial Q^2(p;\lambda ,\theta )}{\partial p^2}& =-\frac{\partial f(x;\lambda ,\theta )/\partial x}{f^2(x;\lambda ,\theta )}\hfill \end{array}.$$

(14)

where $p\ne 1/2$. Further, if we denote:

$$Q_1(p;\lambda ,\theta ):={\left(2p\right)}^{\frac{1}{\theta \lambda }}{\left({\left(2p\right)}^{1/\lambda }+1\right)}^{-\frac{1}{\theta }},Q_2(p;\lambda ,\theta ):={\left({\left(2(1-p)\right)}^{1/\lambda }+1\right)}^{-\frac{1}{\theta }},$$

then, after some computation, it is obtained:

$$\begin{array}{cc}\hfill \frac{\partial Q_1^2(p;\lambda ,\theta )}{\partial p^2}& =-\frac{2^{\frac{1}{\theta \lambda }}{p}^{\frac{1}{\theta \lambda }-2}\left(\theta 2^{1/\lambda }(\lambda +1)p^{1/\lambda }+\theta \lambda -1\right)}{{\theta }^2{\lambda }^2{\left(2^{1/\lambda }{p}^{1/\lambda }+1\right)}^{2+\frac{1}{\theta }}}\hfill \\ \hfill \frac{\partial Q_2^2(p;\lambda ,\theta )}{\partial p^2}& =\frac{{(2-2p)}^{1/\lambda }\left((\theta \lambda +1){(2-2p)}^{1/\lambda }+\theta (\lambda -1)\right)}{{\theta }^2{\lambda }^2{(1-p)}^2{\left({(2-2p)}^{1/\lambda }+1\right)}^{2+\frac{1}{\theta }}}.\hfill \end{array}$$

(15)

Obviously, the sign of derivatives in Equation (15) depends on the sign of the functions:

$$\begin{array}{cc}\hfill {\xi }_1(x;\lambda ,\theta )& =-\theta 2^{1/\lambda }(\lambda +1)p^{1/\lambda }-\theta \lambda +1,\hfill \\ \hfill {\xi }_2(x;\lambda ,\theta )& =(\theta \lambda +1){(2-2p)}^{1/\lambda }+\theta (\lambda -1),\hfill \end{array}$$

which are, respectively, monotonically decreasing and increasing on $p\in (0,1)$. Thus, according to the second one in Equations (14), the RV X has a local maxima at $x_0=2^{-1/\theta }$ if and only if the following inequalities hold:

$${\xi }_1(1/2;\lambda ,\theta )=1-\theta (1+2\lambda )\le 0,{\xi }_2(1/2;\lambda ,\theta )=1-\theta (1-2\lambda )\ge 0.$$

(16)

Therefore, Equations (16) obviously provide a statement of this part of the theorem. □

Remark 4. 

Theorem 4, along with the previously proved Theorem 1, give a complete insight into the variety of shapes of the LLU distribution (see again Figure 1). At the same time, note that its median and mode are equal, but that the true symmetry holds only for $\theta =1$. This can also be confirmed by the definition of the PDF of LLU distribution, given by Equation (2), based on which it follows $f(x;\lambda ,1)=f(1-x;\lambda ,1)$ for each $x\in (0,1)$.

Remark 5. 

Using the QF $Q(p;\lambda ,\theta )$, given by Equation (13), some more measures of shape of the CLU distribution can be studied. These are, for instance, Galton’s skewness (GS) which measures the degree of long tail, and Moors kurtosis (MK) which measures the degree of weight of the tail of the distribution (see, e.g. [12]). In the case of the LLU distribution, these two measures, after some calculation, can be expressed as:

$$\begin{array}{cc}\hfill GS(\lambda ,\theta )& =\frac{Q\left(3/4;\lambda ,\theta \right)-2Q\left(1/2;\lambda ,\theta \right)+Q\left(1/4;\lambda ,\theta \right)}{Q\left(3/4;\lambda ,\theta \right)-Q\left(1/4;\lambda ,\theta \right)}\hfill \\ \hfill & =\frac{2^{-1/\theta }\left(-2^{\frac{1}{\theta \lambda }+1}{\left(2^{-1/\lambda }+1\right)}^{1/\theta }+2^{\frac{\lambda +1}{\theta \lambda }}+2^{1/\theta }\right)}{2^{\frac{1}{\theta \lambda }}-1},\hfill \\ \hfill MK(\lambda ,\theta )& =\frac{Q\left(7/8;\lambda ,\theta \right)-Q\left(5/8;\lambda ,\theta \right)+Q\left(3/8;\lambda ,\theta \right)-Q\left(1/8;\lambda ,\theta \right)}{Q\left(3/4;\lambda ,\theta \right)-Q\left(1/4;\lambda ,\theta \right)}\hfill \\ \hfill & =2^{-\frac{1}{\theta \lambda }}{\left(2^{-1/\lambda }+1\right)}^{1/\theta }\left(\left(2^{\frac{1}{\theta \lambda }}+1\right){\left(4^{-1/\lambda }+1\right)}^{-1/\theta }+\frac{\left(3^{\frac{1}{\theta \lambda }}-4^{\frac{1}{\theta \lambda }}\right){\left({\left(\frac{3}{4}\right)}^{1/\lambda }+1\right)}^{-1/\theta }}{2^{\frac{1}{\theta \lambda }}-1}\right),\hfill \end{array}$$

and 3D plots of their dependence in relation to $(\lambda ,\theta )$ are given in Figure 4.

3. Parameters Estimation & Simulation Study

In this section, the procedure for estimating the unknown parameters $\lambda ,\theta >0$ of the LLU-distributed RV X is presented. In doing so, note that according to the previously described properties of the LLU distribution, using some common procedures to estimate its parameters has some difficulties. Thus, for instance, due to the fact that the moments of the LLU distribution, given by Equation (6), are expressed using the beta function, the method of moments requires certain complex calculations. On the other hand, since the PDF of the LLU distribution is not a differentiable at $x_0=2^{-1/\theta }$, the usage of the maximum likelihood (ML) estimation method is very specific. Considering that the ML estimator of the scale parameter $\lambda$ in the case of the Laplace distribution is equal to the sample median (see, for more detail Norton [26]), here we consider methods of parameter estimation based on the quantiles of the LLU distribution.

For that cause, let $X_1,X_2,\cdots ,X_n$ be the random sample of the length n, for which we define the corresponding order statistics $X_{\left(1\right)}\le X_{\left(2\right)}\le \cdots X_{\left(n\right)}$. As is well-known, PDF of i-th order statistics $X_{\left(i\right)}$ is given as follows:

$$f_{X_{\left(i\right)}}(x;\lambda ,\theta )=\frac{n!}{(i-1)!(n-i)!}f(x;\lambda ,\theta ){\left[F(x;\lambda ,\theta )\right]}^{i-1}{\left[F(x;\lambda ,\theta )\right]}^{n-i},$$

(17)

where $i=1,\cdots ,n$. By substituting $p=j/4$, $j=1,2,3$ into the QF $Q_p:=Q(p;\lambda ,\theta )$, given by Equation (13), for the quartiles of the RV X one obtains:

$$Q(1/4;\lambda ,\theta )=2^{-\frac{1}{\theta \lambda }}{\left(2^{-1/\lambda }+1\right)}^{-1/\theta },Q(3/4;\lambda ,\theta )={\left(2^{-1/\lambda }+1\right)}^{-1/\theta },$$

(18)

while the second quartile is the median $m=Q(1/2;\lambda ,\theta )=2^{-1/\theta }$. On the other hand, the sample quantiles are given by equality:

$${\widehat{Q}}_p:=X_{\left(i_p\right)},i_p=\left\{\begin{array}{cc}np,\hfill & np\mathrm{is}\mathrm{integer},\hfill \\ 1+\left[np\right],\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(19)

where $\left[np\right]$ is the integer part of $np$. In this way, the sample quantiles are in fact the order statistics, and their distribution can be obtained according to Equation (17).

Furthermore, the sample quartiles ${\widehat{Q}}_{j/4}$, $j=1,2,3$ can be used as the estimators of theoretical ones. Thus, by equating median $m=2^{-1/\theta }$ of the RV X with the sample one $\widehat{m}={\widehat{Q}}_{1/2}$, the estimator of the shape parameter $\widehat{\theta }=-ln2/ln\widehat{m}$ is easily obtained. In addition, according to Equations (18) it follows:

$${\left(\frac{Q\left(3/4,\lambda ,\theta \right)}{Q\left(1/4,\lambda ,\theta \right)}\right)}^{\theta }=2^{1/\lambda },$$

and using the estimator $\widehat{\theta }$, for the estimator of the scale parameter $\lambda$ one obtains:

$$\widehat{\lambda }=\frac{ln2}{\widehat{\theta }\left(ln{\widehat{Q}}_{3/4}-ln{\widehat{Q}}_{1/4}\right)}.$$

(20)

In the following, we examine the asymptotic properties of the proposed estimators:

Theorem 5. 

Statistics $(\widehat{\lambda },\widehat{\theta })$ are consistent and asymptotic normal (AN) estimators of the true parameters $(\lambda ,\theta )$, respectively.

Proof. 

First, by using some general asymptotic sample quantile theory [27], we shall prove the consistency of the proposed estimators. For that cause, let us notice that the CDF $p=F(x;\lambda ,\theta )$ is a differentiable and increasing on $p\in (0,1)$. Thus, the quantiles $Q_p:=Q(p;\lambda ,\theta )$ are uniquely determined by Equation (13), and sample quantiles ${\widehat{Q}}_p$ are uniquely determined by Equation (19). Using Bahadur’s representation of sample quantiles (see, e.g. Theorem 1 in Dudek Kuczmaszewska [28], or Serfling [29], pp. 91-92), one obtains:

$${\widehat{Q}}_p=Q_p+\frac{p-F_n\left(Q_p\right)}{f(Q_p;\lambda ,\theta )}+\mathcal{O}\left[{\left(\frac{lnn}{n}\right)}^{3/4}\right],$$

(21)

where $F_n\left(x\right):=n^{-1}{\sum }_{i=1}^{n}1(X_i<n)$ is the empirical CDF. As is known, for each $x\in \mathbb{R}$, the empirical CDF $F_n\left(x\right)$ almost surely and uniformly converge to the CDF $p=F(x;\lambda ,\theta )$, when $n\to \infty$. So, by applying this convergence on Equation (21), when $x=Q_p$, it follows:

$${\widehat{Q}}_p\stackrel{as}{⟶}{Q}_p,n\to \infty ,$$

i.e. the sample quantiles are consistent estimators of the theoretical ones. Finally, as estimators $(\widehat{\lambda },\widehat{\theta })$ are continuous functions of sample quartiles ${\widehat{Q}}_{j/4}$, $j=1,2,3$, by using the property of continuity of almost sure convergence (see, e.g. Serfling [29], pp. 24), it is obtained:

$$(\widehat{\lambda },\widehat{\theta })\stackrel{as}{⟶}(\lambda ,\theta ),n\to \infty .$$

For the proof of the AN property, notice that under the same assumptions sa above, Equation (21) gives the following convergence in distribution:

$$\sqrt{n}\left({\widehat{Q}}_p-Q_p\right)\stackrel{d}{⟶}\mathcal{N}\left(0,\frac{p(1-p)}{{\left[f(Q_p,\lambda ,\theta )\right]}^2}\right),n\to \infty .$$

(22)

According to Equation (22), for the sample median $\widehat{m}={\widehat{Q}}_{1/2}$ one obtains:

$$\sqrt{n}\left(\widehat{m}-m\right)\stackrel{d}{⟶}\mathcal{N}\left(0,{\sigma }_m^2\right),n\to \infty ,$$

wherein is:

$${\sigma }_m^2=\frac{1}{4{\left[f(m;\lambda ,\theta )\right]}^2}{|}_{m=2^{-1/\theta }}={\displaystyle \frac{4^{-\frac{1}{\theta }-1}}{{\lambda }^2{\theta }^2}}.$$

Now, using the continuity of convergence in distribution (see, e.g. Serfling [29], pp. 118), for the estimator $\widehat{\theta }=-ln2/ln\widehat{m}$ is obtained:

$$\sqrt{n}\left(\widehat{\theta }-\theta \right)\stackrel{d}{\to }\mathcal{N}\left(0,{\sigma }_{\theta }^2\right),n\to \infty ,$$

(23)

wherein is:

$${\sigma }_{\theta }^2:=\mathrm{Var}\left(\widehat{\theta }\right)={\left[-\frac{\partial }{\partial m}\left(\frac{ln2}{lnm}\right)\right]}^2{|}_{m=2^{-1/\theta }}·{\sigma }_m^2=\frac{{\theta }^2}{4{\lambda }^2{ln}^{2}2}.$$

The AN property of the estimator $\widehat{\lambda }$, given by Equation (20), is similarly proven. For this purpose, let us first define a statistic:

$$\widehat{\eta }:=\widehat{\theta }\left(ln{\widehat{Q}}_{3/4}-ln{\widehat{Q}}_{1/4}\right),$$

which is consistent estimator of $\eta :=\theta ln\left(Q_{3/4}/Q_{1/4}\right)$. Using Equation (22), as well as the continuity of convergence in distribution, it follows:

$$\sqrt{n}\left(\widehat{\eta }-\eta \right)\stackrel{d}{⟶}\mathcal{N}\left(0,{\sigma }_{\eta }^2\right),n\to \infty ,$$

where, after some computations, one obtains:

$$\begin{array}{cc}\hfill {\sigma }_{\eta }^2& :=\mathrm{Var}\left(\widehat{\eta }\right)={\theta }^2\left(\frac{\mathrm{Var}\left({\widehat{Q}}_{1/4}\right)}{Q_{1/4}^2}+\frac{\mathrm{Var}\left({\widehat{Q}}_{3/4}\right)}{Q_{3/4}^2}\right)=\frac{3{\theta }^2}{16}\left(\frac{1}{Q_{1/4}^2{f}^2\left(Q_{1/4},\lambda ,\theta \right)}+\frac{1}{Q_{3/4}^2{f}^2\left(Q_{3/4},\lambda ,\theta \right)}\right)\hfill \\ \hfill & =\frac{3{\theta }^2}{16}\left(\frac{2^{\frac{2}{\lambda }+4}}{{\theta }^2{\lambda }^2{\left(2^{1/\lambda }+1\right)}^2}+\frac{16}{{\theta }^2{\lambda }^2{\left(2^{1/\lambda }+1\right)}^2}\right)=\frac{3\left(4^{1/\lambda }+1\right)}{{\lambda }^2{\left(2^{1/\lambda }+1\right)}^2}.\hfill \end{array}$$

(24)

Finally, according to $\widehat{\lambda }=ln2/\widehat{\eta }$ and Equation (24), it is obtained:

$$\sqrt{n}\left(\widehat{\lambda }-\lambda \right)\stackrel{d}{⟶}\mathcal{N}\left(0,{\sigma }_{\lambda }^2\right),n\to \infty ,$$

(25)

wherein is:

$${\sigma }_{\lambda }^2:=\mathrm{Var}\left(\widehat{\lambda }\right)={\left[\frac{\partial }{\partial \eta }\left(\frac{ln2}{\eta }\right)\right]}^2·{\sigma }_{\eta }^2=\frac{3{\lambda }^2\left(4^{1/\lambda }+1\right)}{{ln}^{2}2{\left(2^{1/\lambda }+1\right)}^2}.$$

(26)

Thus proven convergences in Equations (23) and (25) confirm the AN properties of both proposed estimators $(\widehat{\theta },\widehat{\lambda })$. □

Remark 6. 

It is worth noting that the variance of the estimator $\widehat{\lambda }$, given by Equation (26), does not depend on the parameter θ. This, among others, is one of the reasons that justifies the use of this estimator.

A numerical study of the effectiveness of the proposed estimators is presented below, based on Monte Carlo simulations of samples $x_1,\cdots ,x_n$ taken from the LLU distribution. More precisely, for different samples and parameters values, the proposed estimators were calculated and their statistical analysis was performed. To this end, three different types of samples are considered (see also Figure 5 below):

$\left(i\right)$ Sample I is drawn from an increasing LLU distribution with parameters $\lambda =0.3$ and $\theta =3.5$.

$\left(ii\right)$ Sample II was drawn from a symmetric, both-sides tailed (i.e. bathtub-shaped) LLU distribution with parameters $\lambda =0.3$ and $\theta =1$.

$\left(iii\right)$ Sample III was drawn from a symmetric, both-sides vanishing, unimodal LLU distribution with parameters $\lambda =2$ and $\theta =1$.

The simulated values of all samples are generated by the R-package "distr" [30], according to which the estimates $\widehat{\theta }$ and $\widehat{\lambda }$ are calculated. In order to additionally check the effectiveness of the proposed estimators, the sample realizations of different lengths $n\in \{150,500,1500\}$ were observed. It is worth noting that the above sample sizes were chosen similar to some of the real-world data which will be further analyzed. Additionally, $S=200$ independent simulations were performed for each sample, and the results of their statistical analysis are presented in the following Table 1, Table 2 and Table 3.

More precisely, Table 1, Table 2 and Table 3 contain the summary statistics of the obtained estimates, that is, their minimums (Min.), mean values (Mean), maximums (Max.) and standard deviations (SD). Additionally, mean square estimation errors (MSEE), fractional estimation errors (FEE), and the Anderson-Darling (AD) normality test statistics are provided, too. According to this, it can be observed that the proposed estimators are efficient, as the bias, sample range (Max.–Min.), as well as SD and MSEE values decrease with increasing sample size. Thereby, it can be noticed that the estimates $\widehat{\theta }$ are more stable and efficient than $\widehat{\lambda }$. This is due to the fact that the estimate $\widehat{\lambda }$ is obtained by a two-step procedure, according to the previously calculated values of $\widehat{\theta }$, and then using Equation (20).

Similar conclusions can be confirmed based on the AN analysis of the obtained estimates. Its examination was performed using the Anderson-Darling normality test, whose test statistic (AD) and corresponding p-values are calculated using the R-package "nortest" [31]. According to thus obtained results, also presented in Table 1, Table 2 and Table 3, it can be noted that estimates of $\widehat{\lambda }$ from smaller samples have a less pronounced AN feature than $\widehat{\theta }$. Nevertheless, the AN property was confirmed in most cases, especially for larger samples. Some confirmation of these facts can be seen in Figure 5, where realizations, empirical and theoretical distributions of the observed samples are shown.

4. Applications of the LLU distribution

This section discusses some possibilities of applying the LLU distribution in real-world data modeling, primarily in the domain of dynamic and regression analysis. Also, a comparison of the LLU distribution with some existing, well-known and frequently used unit distributions is made. For this purpose, three sets of data are considered, and their brief description is as follows:

$\left(i\right)$ The first data set, called Series A, represents broadband usage in $n=153$ rural counties in the United States, based on Microsoft’s Air Belt initiative to help close the rural broadband gap and improve the performance and security of broadband software and services. More specifically, this dataset, taken from the GitHub, Inc. database [32], consists of the percentage of devices connected to the Internet at broadband speed by each zip code, during October 2020.

$\left(ii\right)$ The second one (Series B), taken from the official website of the National Center for Environmental Information, contains historical data on the melting rate of the South Greenland ice core. The dataset itself was based on research and reconstruction of ice melting rates conducted by Kameda et al. [33] and the World Data Center for Paleoclimatology. In this way, a time series of annual percentage data was generated in the period from 1546. to 1989., the length of which is $n=444$.

$\left(iii\right)$ Finally, the third, Series C, is obtained according to the official data from the National Association of Securities Dealers Automated Quotations (NASDAQ) Stock Market [34], and the so-called log-returns of daily changes of natural gas prices (in US dollars per cubic meter), from 1 January, 2018. until to 1 March, 2023. It represents a time series of length $n=1302$, which was also considered in Stojanović et al. [35], where it was shown that it can be viewed as a series of independent realizations of RVs with Laplace distribution. Therefore, this series was transformed using the previously mentioned logistic function $x=\phi (y;\theta )$ and thus is obtained a corresponding regression model with an output variable that can be seen as the realization of independent RVs with LLU distribution.

In order to additionally verify the efficiency of the LLU distribution, it is compared with the well-known Beta distribution and the Kumaraswamy distribution, whose PDFs are, respectively,

$$g_1(x;a,b)=\frac{1}{B(a,b)}{x}^{a-1}{(1-x)}^{b-1},g_2(x;a,b)=ab{x}^{a-1}{(1-x^a)}^{b-1}.$$

Herein is $0<x<1$, while $a,b>0$ are the distribution parameters, and $B(a,b)$ is the beta function. To obtain the estimated values of the parameters of the Beta distribution, the method of moments (MM) is applied, according to which the estimates are as follows:

$${\displaystyle \widehat{a}=\overline{x}\left(\frac{\overline{x}(1-\overline{x})}{\overline{v}}-1\right),\widehat{b}=(1-\overline{x})\left(\frac{\overline{x}(1-\overline{x})}{\overline{v}}-1\right).}$$

Here, $\overline{x}=n^{-1}{\sum }_{i=1}^n{x}_i$ and $\overline{v}={(n-1)}^{-1}{\sum }_{i=1}^n{(x_i-\overline{x})}^2$ are the sample mean and variance, respectively, wherein is $\overline{v}<\overline{x}(1-\overline{x})$. On the other hand, for the Kumaraswamy distribution, the maximum likelihood (ML) estimation method is used, according to which the following parameter estimates are simply obtained (see, e.g. Jones [36]):

$$\widehat{a}=-{\left(\frac{1}{n}\sum _{i=1}^{n}ln{x}_i\right)}^{-1},\widehat{b}=-{\left(\frac{1}{n}\sum _{i=1}^{n}ln\left(1-x_i^{\widehat{a}}\right)\right)}^{-1}.$$

We emphasize that one of the reasons for the choice of estimation methods of these two distributions is the comparison not only with their distributions, but also with regard to different estimation procedures.

Realizations of the series mentioned above, along with the fitted PDFs obtained by the previously described estimation procedures, are shown in Figure 6. As can be seen, all series have pronounced and persistent fluctuations, which create "heavy tails" in their distributions. Specifically, Series A has an increasing PDF, for Series B the PDF is positively asymmetric and unimodal, while the PDF of Series C indicates the symmetry of its distribution. Estimated parameter values for each series and for all competing models are shown in Table 4, and $S=1000$ independent Monte Carlo simulations were conducted using them. Thereafter, the agreement of the distributions between the actual and fitted data was checked using the MSEE error statistic, the Akaike information criterion (AIC ), and the Kolmogorov- Smirnov (KS) test of equivalence of the asymptotic distribution of two samples. Based on the results thus obtained, it is noticeable that MSEE and AIC values are generally lower when LLU and Beta distribution are applied as the appropriate fitting model. At the same time, it is obvious that the LLU distribution has better fit-characteristics than both other distributions. Moreover, only in the case of the LLU distribution the KS test do not reject, with a significant level $p>0.01$, the hypothesis of equivalence with the observed empirical distributions.

5. Conclusion

A new, the so-called the LLU distribution is presented here, along with the corresponding key properties and a procedure for estimating its parameters based on quantiles. The consistency and AN property of the estimators were also examined, as well as a Monte-Carlo study of their efficiency. Finally, a practical application of the LLU distribution in fitting real-world data is presented. To verify the effectiveness of the proposed model, it was applied to fit the distributions of three real-world datasets and compared with the Beta and Kumaraswamy distributions. According to the results obtained in this way, the LLU distribution provides a better fit to the observed data, which is a motivation for further research of some other new unit distributions, with slightly different characteristics. Thus, for instance, by applying the logistic map to some other probability distribution supported on the entire real line (such as, for example, normal, Student or Gumbel distribution), new unit distributions can be obtained and this can be a guideline for some potential further researches for the authors.