A New Class of Probability Distributions via Bi-Elliptic Functions

1. Introduction

The topic of continuous distributions on a bounded interval is an important aspect of statistical theory and applications. In recent years, many researchers have focused their attention on developing such distributions and their applications (Kotz and Dorp, 2004; García et al., 2011; Condino and Domma, 2017). Beta distribution is one of the most popular bounded distributions and is useful in many areas of application because it provides a rich family of distributional shapes over some finite interval (Goel and Strebel, 1984; McDonald and Xu, 1995; Gupta and Nadarajah, 2004). Researchers developed suitable substitutes of the beta distribution over the years. For instance, the triangular distribution has been investigated by D. Johnson (1997) as a proxy for the beta distribution in risk analysis. Many extensions have been proposed over the years. The trapezoidal distribution was constructed by adding an additional horizontal line segment in the middle and thus provides a more flexible distribution. It was used to model the duration and the form of a phenomenon which maybe represented by three stages: a growth-stage, a stable stage, and a decline stage in Dorp and Kotz (2003). Dorp and Kotz (2002) developed a new family of distributions named the two-sided power distribution and discussed its parameter estimation. Karlis and Xekalaki (2008) introduced another class of distributions called the polygonal distributions constructed from finite mixtures of the triangular distributions. Vander Wielen and Vander Wielen (2015) proposed the general segmented distributions constructed from a finite number of line segments.

In this paper, a new class of distributions is developed using two elliptical arcs connecting together over a bounded interval. We term this distribution the bi-elliptic distribution. Similar to the triangular distribution, it is a unimodal distribution with three parameters: the minimal value a, the maximal value b, and the mode m. The proposed distribution has several advantages in probability modeling. First, it is differentiable at all points over its support. Unlike the triangular distribution and other similar truncated distributions, the density function of the bi-elliptic distribution is differential everywhere in its domain, making it a better choice analytically. Second, the proposed distribution has a flatter shape around its mode, making it a better modeling choice for certain physical and social science phenomena of interest. For instance, small businesses usually experience slow growth and decline when they reach the stage of maturity (Scott and Bruce,1987). Third, the three parameters in the distribution provide large flexibility in modeling. For instance, the distribution can be negatively skewed, positively skewed, and symmetric. Finally, similar to the triangular distribution, the bi-elliptic distribution is based on a knowledge of the minimum, maximum, and the most common modal value, enabling it to work when data are limited. For instance, it can be used in quantitatively analyzing the uncertainty of three-point-estimate problems (Mulligan, 2016).

The paper is organized as follows. Section 2 presents the bi-elliptic distribution and its basic properties. In Section 3, the maximum likelihood estimation (MLE) procedure for bi-elliptic distributions is discussed. Section 4 presents a wide range of Monte-Carlo simulation results for parameter estimation. Section 5 compares the triangular distribution and the proposed bi-elliptic distribution in fitting data drawn from a range of beta distributions. Finally, we provide some concluding remarks in Section 6.

2. The Bi-Elliptic Distribution

Let X be a random variable with a probability density function given by

$$f(x;a,b,m)=\left\{\begin{array}{cc}{\displaystyle \frac{h}{r_1}}\sqrt{r_1^2-{(x-m)}^2},\hfill & a\le x\le m\hfill \\ {\displaystyle \frac{h}{r_2}}\sqrt{r_2^2-{(x-m)}^2},\hfill & m<x\le b\hfill \end{array}\right.$$

(1)

where $r_1=m-a$, $r_2=b-m$, and $h={\displaystyle \frac{4}{\pi (b-a)}}$. The random variable X is said to have a bi-elliptic distribution, $E(a,b,m)$, $a<m<b$. Here a is the lower limit, b is the the upper limit, and m is the mode. In addition, bi-elliptic distributions also include two extreme cases when $m=a$ or $m=b$:

$$\begin{array}{cc}\hfill f(x;a,b,m)& ={\displaystyle \frac{h}{r_1}}\sqrt{r_1^2-{(x-m)}^2},a\le x\le b,ifm=b,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill f(x;a,b,m)& ={\displaystyle \frac{h}{r_2}}\sqrt{r_2^2-{(x-m)}^2},a\le x\le b,ifm=a.\hfill \end{array}$$

(3)

We call the bi-elliptic distribution with $a=0$ and $b=1$ the standard bi-elliptic distribution. For any bi-elliptic distribution $X\sim E(a,b,m)$, the following property holds:

$$Y={\displaystyle \frac{X-a}{b-a}}\sim E\left(0,1,{\displaystyle \frac{m-a}{b-a}}\right).$$

(4)

Figure 1 provides the density plots of examples of symmetric, positively and negatively skewed $E(0,1,m)$ distributions for different choices of m values. As we can see, the distribution is right-skewed if $m<0.5$; it is left-skewed if $m>0.5$; the distribution is symmetric if $m=0.5$.

The cumulative distribution function of an $E(a,b,m)$ distribution can be derived from equation (1) and is given below

$$F\left(x\right)=\left\{\begin{array}{cc}{\displaystyle \frac{h}{2r_1}}\left[(x-m)\sqrt{r_1^2-{(x-m)}^2}+r_1^2{sin}^{-1}\left({\displaystyle \frac{x-m}{r_1}}\right)\right]+{\displaystyle \frac{\pi r_{1}h}{4}},\hfill & a\le x\le m\hfill \\ {\displaystyle \frac{h}{2r_2}}\left[(x-m)\sqrt{r_2^2-{(x-m)}^2}+r_2^2{sin}^{-1}\left({\displaystyle \frac{x-m}{r_2}}\right)\right]+{\displaystyle \frac{\pi r_{1}h}{4}},\hfill & m<x\le b.\hfill \end{array}\right.$$

(5)

Figure 2 provides plots of the distribution function for examples of symmetric, positively and negatively skewed $E(0,1,m)$ distributions for different choices of m values.

For any bi-elliptic distribution, m is the mode and the $\alpha$-quantile $Q_{\alpha }\left(x\right)$ with $\alpha ={\displaystyle \frac{m-a}{b-a}}$.

The median of a bi-elliptic distribution $E(a,b,m)$ is the unique real value x for which $F\left(x\right)=0.5$. There is no general closed-form expression for the median. However, for symmetric cases when $m=(a+b)/2$, the median is the mode m.

The first three moments and the variance of $E(a,b,m)$ can be obtained from the density function (1) and simplified as follows:

$$\begin{array}{cc}\hfill E\left(X\right)& =\mu =m+{\displaystyle \frac{4(a+b-2m)}{3\pi }},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill E\left(X^2\right)& =m^2+{\displaystyle \frac{8m(a+b-2m)}{3\pi }}+{\displaystyle \frac{{(m-a)}^2-(m-a)(b-m)+{(b-m)}^2}{4}},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill E\left(X^3\right)& =m^3+{\displaystyle \frac{4}{\pi (b-a)}}\left[{\displaystyle \frac{2}{15}}(r_2^4-r_1^4)+m^2(r_2^2-r_1^2)\right]\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & +{\displaystyle \frac{4}{(b-a)}}\left[{\displaystyle \frac{3m}{16}}(r_1^3+r_2^3)+{\displaystyle \frac{m^3}{4}}(r_1+r_2)\right]\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill Var\left(X\right)& ={\sigma }^2={\displaystyle \frac{{(m-a)}^2-(m-a)(b-m)+{(b-m)}^2}{4}}-{\displaystyle \frac{16{(a+b-2m)}^2}{9{\pi }^2}}.\hfill \end{array}$$

(10)

Note that the skewness of X can be expressed in term of the first three moments of X:

$$\begin{array}{c}\hfill skew\left(X\right)={\displaystyle \frac{E\left(X^3\right)-3\mu {\sigma }^2-{\mu }^3}{{\sigma }^3}}.\end{array}$$

(11)

3. Parameter Estimation

3.1. Maximum Likelihood Estimation

For a random sample $\mathbf{X}=(X_1,\cdots ,X_n)$ of size n from an $E(a,b,m)$ distribution, let the order statistics be $X_{\left(1\right)}<X_{\left(2\right)}<\cdots <X_{\left(n\right)}$. Using expression (1), the likelihood function for $\mathbf{X}$ can be obtained as follows:

$$L(\mathbf{X};a,b,m)={\left[{\displaystyle \frac{4}{\pi (b-a)}}\right]}^{n}H(\mathbf{X};a,b,m),$$

(12)

where

$$H(\mathbf{X};a,b,m)=\prod _{i=1}^s\sqrt{1-{\left({\displaystyle \frac{X_{\left(i\right)}-m}{m-a}}\right)}^2}·\prod _{i=s+1}^n\sqrt{1-{\left({\displaystyle \frac{X_{\left(i\right)}-m}{b-m}}\right)}^2},$$

(13)

$X_{\left(0\right)}\equiv a,X_{(n+1)}\equiv b$ and s is implicitly defined by $X_{\left(s\right)}\le m\le X_{(s+1)}.$

From equation (12), the log-likelihood function can be written as

$$logL(\mathbf{X};a,b,m)=nlog\left({\displaystyle \frac{4}{\pi (b-a)}}\right)+{\displaystyle \frac{1}{2}}\sum _{i=1}^{s}log\left(1-{\left({\displaystyle \frac{X_{\left(i\right)}-m}{m-a}}\right)}^2\right)+{\displaystyle \frac{1}{2}}\sum _{i=s+1}^{n}log\left(1-{\left({\displaystyle \frac{X_{\left(i\right)}-m}{b-m}}\right)}^2\right),$$

(14)

which is equivalent to:

$$\begin{array}{cc}\hfill logL(\mathbf{X};a,b,m)& =nlog{\displaystyle \frac{4}{\pi }}-nlog(b-a)\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\sum _{i=1}^s\left[log\left(X_{\left(i\right)}-a\right)+log\left(2m-a-X_{\left(i\right)}\right)-2log(m-a)\right]\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\sum _{i=s+1}^n\left[log\left(b-X_{\left(i\right)}\right)+log\left(b+X_{\left(i\right)}-2m\right)-2log(b-m)\right].\hfill \end{array}$$

(15)

We have the following partial derivatives:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial logL}{\partial a}}& ={\displaystyle \frac{n}{b-a}}+{\displaystyle \frac{1}{2}}\sum _{i=1}^s\left[{\displaystyle \frac{1}{a-X_{\left(i\right)}}}+{\displaystyle \frac{1}{a+X_{\left(i\right)}-2m}}-{\displaystyle \frac{2}{a-m}}\right],\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial logL}{\partial b}}& =-{\displaystyle \frac{n}{b-a}}+{\displaystyle \frac{1}{2}}\sum _{i=s+1}^n\left[{\displaystyle \frac{1}{b-X_{\left(i\right)}}}+{\displaystyle \frac{1}{b+X_{\left(i\right)}-2m}}-{\displaystyle \frac{2}{b-m}}\right],\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial logL}{\partial m}}& =-\sum _{i=1}^s\left[{\displaystyle \frac{1}{a+X_{\left(i\right)}-2m}}-{\displaystyle \frac{1}{a-m}}\right]-\sum _{i=s+1}^n\left[{\displaystyle \frac{1}{b+X_{\left(i\right)}-2m}}-{\displaystyle \frac{1}{b-m}}\right].\hfill \end{array}$$

(18)

3.2. One-Parameter Maximum Likelihood Estimation

The one-parameter case is of special interest because in some applications, there exist natural lower and upper bounds for random variables of interests. For instance, in hydrology, the reservoir yield and storage distribution has an natural upper bound being its capacity and a lower bound being zero (Fletcher and Ponnambalam, 1996).

When $a=a_0$ and $b=b_0$ are fixed and known, taking the partial derivative of the log-likelihood function $logL$ with respect to m yields equation (18), whose critical point(s) do not have an explicit form. The one-parameter maximum likelihood estimation is an optimization problem below:

$$\underset{a_0<m<b_0}{max}\left\{L(\mathbf{X};a_0,b_0,m)\right\}$$

(19)

or, equivalently,

$$\underset{a_0<m<b_0}{max}\left\{logL(\mathbf{X};a_0,b_0,m)\right\}$$

(20)

where $L(\mathbf{X};a_0,b_0,m)$ and $logL(\mathbf{X};a_0,b_0,m)$ are given by equations (12) and (14), respectively, with a and b being replaced by $a_0$ and $b_0$. This optimization problem can be solved using numerical methods. For instance, we use Brent’s algorithm for the numerical computations which combines linear interpolation and inverse quadratic interpolation with bisection for one-dimensional optimization problems. Interested readers may refer to Brent (1971) for more details.

3.3. Three-Parameter Maximum Likelihood Estimation

There are many cases where we have no information about the values of a, b, and m, and in such cases we need to estimate the values of the three parameters simultaneously. The partial derivatives of the log likelihood function $logL(\mathbf{X};a,b,m)$ are derived and given in equations (16), (12), and (14), respectively. Again, we cannot find explicit forms for their critical values and we resort to numerical methods for finding the MLE via the log likelihood function. We have

$$\underset{S(a,b,m)}{max}\left\{logL(\mathbf{X};a,b,m)\right\}$$

(21)

where $S(a,b,m)=\{(a,m,b):a<X_{\left(1\right)},b>X_{\left(n\right)},a<m<b\}$ and $logL(\mathbf{X};a,b,m)$ is given by equation (14). A linearly constrained optimization algorithm can be used for the numerical computation. Given a random sample, we can calculate gradients (vector of the partial derivatives with respect to a, b, and m) of the log likelihood function $logL(\mathbf{X};a,b,m)$. For a random sample $\mathbf{X}=(X_1,\cdots ,X_n)$ with order statistics $X_{\left(1\right)}<X_{\left(2\right)}<\cdots <X_{\left(n\right)}$. The linear constraints on the parameters a, b , and m can be formulated as follows

$$a<X_{\left(1\right)},b>X_{\left(n\right)},m>a,m<b.$$

Interested readers may refer to Lange (2016) for details of the linearly constrained optimization algorithm.

4. Monte Carlo Simulations

In this section we illustrate some previous results on one-parameter MLE and three-parameter MLE using Monte Carlo simulations. All computations are performed using the open source statistical software R.

4.1. One-Parameter MLE

For one-parameter MLE, we consider the bi-elliptic distributions where m is the only unknown parameter. For our simulation, we set $a=0$ and $b=1$.

The values of m and n used in the simulation are $m=0.3,0.5,0.7$; $n=200,400,600,800,1000$. This gives a total of 15 different combinations of m and n. For each fixed pair of m and n, 1000 random samples were drawn from the corresponding bi-elliptic distribution $f(x;0,1,m)$. For each of the 1000 random samples, we calculate the estimates of the maximum likelihood estimator using the numerical method described in the previous section. The initial value is set at zero. These estimates are denoted by ${\widehat{m}}_1,\cdots ,{\widehat{m}}_{1000}$. Four performance measures are used to evaluate how our method estimates the value of the parameter m: the bias, the empirical standard error (Empirical SE), the mean square error (MSE), and the root mean square error (RMSE). Table 1 summarizes these measures.

Table 2 presents the simulation results for the one-parameter MLE. We have two conclusions. First, for each selected value of m, the absolute values of the biases associated with the estimator $\widehat{m}$ in estimating m have a decreasing trend as the sample size increases from 200 to 1000. The Empirical SE, MSE, and RMSE for $\widehat{m}$ in estimating m are also given and they all decrease as the sample size increases. In addition, for each selected sample size n, the values of the four performance measures for different values of m are, in general, close to each other. This implies that our estimator $\widehat{m}$ has a uniform effect in estimating m for different choices of m.

4.2. Three-Parameter MLE

In the 3-parameter MLE, all the three parameters a, b, and m in the bi-elliptic distribution are unknown. For our simulation, we set $a=-1$, $b=5$, and $m=0.3$. The sample sizes used are $n=200,400,600,800,1000$. 1000 random samples were drawn from the bi-elliptic distribution $f(x;-1,5,0.3)$. For each sample, we calculate the estimates $\widehat{a}$, $\widehat{b}$, and $\widehat{m}$ numerically. The initial values are set at: $a_0=a-(b-a)/4=-2.5,b_0=b+(b-a)/4=6.5,$ and $m_0$ is one random value between -2.5 and 6.5 drawn from the uniform distribution $U(-2.5,6.5)$.

These estimates are denoted by ${\{{\widehat{a}}_i,{\widehat{b}}_i,{\widehat{m}}_i\}}_{i=1}^{1000}$. Similar to the 1-parameter MLE, we use the same four performance measures listed in Table 1.

Table 3 presents the simulation results for the three-parameter MLE. The absolute values of the biases have a decreasing trend as the sample size increases. The Empirical SE, MSE, and RMSE for $\widehat{m}$ in estimating m are also given and they all decrease as the sample size increases as well. The values of the performance measures for different samples are close for $\widehat{a}$ and $\widehat{b}$. However, the values for these four performance measure for $\widehat{m}$ are in general significantly larger than those for $\widehat{a}$ and $\widehat{b}$.

5. Comparing the Bi-Elliptic Distribution and the Triangular Distribution as a Proxy for the Beta Distribution

In the literature, one application of the triangular distribution is to serve as a proxy for the beta distribution. Johnson (1997) noted that “the differences between the two distribution functions are seldom significant". In this section, we present empirical results comparing the bi-elliptic distribution and the triangular distribution as a proxy for the beta distribution. We consider 81 beta distributions resulting from combinations of 9 different values of the two shape parameters $\alpha$ and $\beta$. For each selected beta distribution, we draw 1000 samples of size 100, fit the data with the bi-elliptic distribution and the triangular distribution. Based on the MLE estimates, we compute the corresponding AICs and the Akaike weights of the two models. To obtain the Akaike weights, we first compute, for each model, the differences in AIC with respect to the AIC of the best candidate model (Burnham and Anderson, 2002):

$${\Delta }_i\left(AIC\right)=AI{C}_i-min\left(AIC\right),$$

where i is either “bi-elliptic" or “triangular".

Then the Akaike weights are computed:

$$w_i\left(AIC\right)={\displaystyle \frac{exp\{-\frac{1}{2}{\Delta }_i\left(AIC\right)\}}{\sum exp\{-\frac{1}{2}{\Delta }_k\left(AIC\right)\}}}.$$

Weight $w_i\left(AIC\right)$ can be interpreted as the probability that $M_i$ is the best model (in the AIC sense, that it minimizes the Kullback–Leibler discrepancy), given the data and the set of candidate models (Burnham and Anderson, 2002).

To compare the bi-elliptic model and the triangular model, we compute the ratio

$$r={\displaystyle \frac{w_{bi-elliptic}}{w_{bi-elliptic}+w_{triangular}}},$$

which gives a the normalized probability that the bi-elliptic model is to be preferred over the triangular model.

Table 4 reports the mean AIC, the percentage of times that AIC value corresponding to the bi-elliptic distribution ($AI{C}_{be}$)is smaller than that of the triangular distribution ($AI{C}_t$) in 1000 simulations, and the ratio r for different combinations of the shape parameters $\alpha$ and $\beta$. Compared to triangular distributions, bi-elliptic distributions offer a significantly better proxy for the beta distributions with different combinations of $\alpha$ and $\beta$ between 1 and 2. For instance, when $\alpha =1$, there is a probability between 81% and 99% that the bi-elliptic proxy is superior compared to the triangular proxy for $\beta$ between 1 and 1.4. As $\alpha$ increases, the range for “acceptable" $\beta$ (in terms of providing a high probability r of at least 70%) broadens, but never exceeds $\beta =2$. Due to the mirror-image symmetry of the beta densities, similar behavior of the probability r can be observed when $\beta$ is fixed and $\alpha$ varies. Figure 3a shows the density of the beta distribution with $\alpha =1.4$ and $\beta =1.2$ with the fitted curves using the bi-elliptical distribution and the triangular distribution. When the beta density resembles the bi-elliptical distribution, the later proves much better proxy compared to the triangular distribution.

On the other hand, when either $\alpha$ or $\beta$ is greater than 2, the triangular distribution appears to be a better proxy for the beta distribution. For example, when $\beta =1$, there is a probability between 0.5% and 20% that the bi-elliptic proxy is superior compared to the triangular proxy for $\alpha$ between 2 and 3.5. Or equivalently, there is a probability between 80% and 99.5% that the triangular proxy is more preferable to the bi-elliptic one for $\alpha$ between 2 and 3.5 and $\beta =1$. As seen in Figure 3b, where the beta density with $\alpha =3$ and $\beta =3.5$ plotted together with the bi-elliptical and triangular fits, triangular fit appears to be a better fit.

6. Conclusions

A new family of bi-elliptic distributions supported on a bounded interval has been introduced by connecting two elliptical arcs. Similar to the triangular distribution, the proposed distribution has three parameters and can take a variety of shapes. These elliptical shapes offer an appealing potential as modeling tools in areas such as risk analysis.

Important properties of the bi-elliptic distribution has been discussed. The MLE algorithms have been introduced for both the one-parameter estimation and the three-parameter estimation problems. A Monte Carlo simulation study has shown satisfying estimation results.

A comparison study was conducted to explore the use of the bi-elliptic distribution versus the triangular distribution as an approximation to the beta distribution. Results show that, for a certain range of beta distributions’ parameter values, the proposed distribution is superior compared to the triangular distribution. The bi-elliptic distribution seems to be a useful competitor to the triangular distribution and the beta distribution, especially in applications such as risk analysis. It is our hope that the introduction of the bi-elliptic distribution into statistical theory and practice will contribute to the basic goals of applied statistical studies.