On the Arcsecant Hyperbolic Normal Distribution. Properties, Quantile Regression Modeling and Applications

Abstract

This work proposes a new distribution defined on the unit interval. It is obtained by a novel transformation of a normal random variable involving the hyperbolic secant function and its inverse. The use of such a function in distribution theory has not received much attention in the literature, and may be of interest for theoretical and practical purposes. Basic statistical properties of the newly defined distribution are derived, including moments, skewness, kurtosis and order statistics. For the related model, the parametric estimation is examined through different methods. We assess the performance of the obtained estimates by two complementary simulation studies. Also, the quantile regression model based on the proposed distribution is introduced. Applications to three real datasets show that the proposed models are quite competitive in comparison to well-established models.

1. Introduction

Over the past twenty years, many statisticians and researchers have focused on proposing new extended or generalized distributions by adding additional parameters to the basic probability distributions. The common point of these studies is to obtain better inferences than those of the baseline probability distributions. In this context, especially, the modeling approaches on the unit interval have recently multiplied since they are related to specific issues such as the recovery rate, mortality rate, daily patient rate, etc. The beta distribution is the best-known distribution defined over the unit interval for modeling the above measures. It has great flexibility in the shapes of the probability density function (pdf) and hazard rate function (hrf). Although it has very flexible forms for data modeling, sometimes it is not sufficient for modeling and explaining unit datasets. For this reason, new alternative unit models have been proposed in the statistical distribution literature, including the Johnson_SB [1], Topp-Leone [2], Kumaraswamy [3], standart two-sided power [4], log-Lindley [5], log-xgamma [6], unit Birnbaum-Saunders [7], unit Weibull [8], unit Lindley [9], unit inverse Gaussian [10], unit Gompertz [11], second degree unit Lindley [12], log-weighted exponential [13], logit slash [14], unit generalized half normal [15], unit Johnson_SU [16], trapezoidal beta [17] and unit Rayleigh [18] distributions. Many of the above distributions were obtained by transforming the baseline distribution, and they performed better than the beta distribution in terms of data modeling. For instance, the Johnson_SBdistribution was created via logistic transformation of the ordinary normal distribution. In this way, a very flexible unit normal distribution was obtained over the unit interval. The other mentioned unit distributions introduced over the last decade can also be seen as alternatives to the well-known beta, Johnson_SB, Topp-Leone and Kumaraswamy distributions.

On the other hand, the ordinary regression models explain the response variable for given certain values of the covariates based on the conditional mean. However, the mean may be affected by a skewed distribution or outliers in the measurements. Possible solutions are provided by the quantile regression models proposed by [19], particularly popular for being less sensitive to outliers than the ordinary regression models.

In line with above, the aim of this study is to introduce a new alternative unit probability distribution based on the normal distribution. More precisely, we use a new transformation of the normal distribution based on the hyperbolic secant function. As a matter of fact, the use of the hyperbolic function has not received enough attention in the published literature on distribution theory, despite the great interest among students and practitioners of the few distributions based on it. Examples include the famous hyperbolic secant distribution and its generalizations as presented in [20]. In a sense, we show that the proposed methodology allows us to transport the applicability and working capacity of the normal distribution to the unit interval. In particular, we develop a new quantile regression modeling via the re-parameterizing of the new probability distribution in terms of any quantile. All these aspects are developed in the article through mathematical, graphical and numerical approaches.

The paper has been set as follows. We define the proposed distribution in Section 2. Its basic distributional properties are described in Section 3. Section 4 is devoted to the procedures of the different parametric estimation methods. Two different simulation studies are given to see the performance of the different estimates of the model parameters in Section 5. The new quantile regression model based on the proposed distribution and its residual analysis are introduced by Section 6. Three real data illustrations, one of which relates to quantile modeling and others to univariate data modeling, are illustrated in Section 7. Finally, the paper is ended with conclusions in Section 8.

2. The New Unit Distribution and Its Properties

The new unit distribution is defined as follows: Let Y be a random variable such thatY∼N(μ,^σ2)whereμ∈Randσ>0, and X be the random variable defined by

X=sechY,

wheresechy=2/(^ey+^e−y)=2^ey/(^e2y+1)∈(0,1)is the hyperbolic secant function fory∈R, also known as the inverse of the hyperbolic cosine function. Then the distribution of X is called “arcsech” normal distribution and it is denoted byASHNorASHN(μ,σ)whenμandσare required. To our knowledge, it constitutes a new unit distribution; It is unlisted in the literature. Before stating the motivations for theASHNdistribution, the corresponding cumulative distribution function (cdf) and pdf are presented in the following proposition.

Proposition 1.

The cdf and pdf of theASHN(μ,σ)distribution are given as

F(x,μ,σ)=2−Φarcsechx+μσ−Φarcsechx−μσ

and

f(x,μ,σ)=1σx1−^x2ϕarcsechx+μσ+ϕarcsechx−μσ,

respectively, forx∈(0,1), wherearcsechz=log1+1−^z2/z>0is the hyperbolic arcsecant function (or inverse hyperbolic secant function) forz∈(0,1),Φ(x)andϕ(x)are the cdf and pdf of theN(0,1)distribution, respectively. Forx∉(0,1), standard completions on these functions are performed.

For the sake of presentation, the proof of this result and those of the results to come are given in Appendix A.

Based on Proposition 1, as a first property, note that, forμ=0andx∈(0,1), the cdf and pdf are reduced to the quite manageable functions:

F(x,0,σ)=21−Φ1σarcsechx

and

f(x,0,σ)=2σx1−^x2ϕ1σarcsechx.

In full generality, forx∈(0,1), an alternative formulation for the pdf is

f(x,μ,σ)=2π1σx1−^{x2e−arcsechx2+μ22σ2}coshμ^σ2arcsechx,

wherecoshy=(^ey+^e−y)/2is the hyperbolic cosine function fory∈R . Eventually, we can express the cosh term in Equation (5) as

coshμ^σ2arcsechx=12^{1x+1x2−1μσ2}+^{1x+1x2−1−μσ2}.

Let us now focus on the behavior off(x,μ,σ)at the boundaries.

• When x tends to 0, sincearcsechx∼−logx→+∞and it appears in power 2 the exponential term, we havef(x,μ,σ)→0.

• When x tends to 1, sincearcsech1=0, we have

  f(x,μ,σ)∼1π1σ1−x^e−μ22σ2→+∞.

  Ifσis large and^μ2≈2^σ2, or^μ2/2^σ2is large, the pointx=1appears as a “special singularity” in the following sense: The functionf(x,μ,σ)can decrease to 0 in the neighborhood ofx=1, then suddenly explodes atx=1. This phenomenon is only punctual; this is not a particular disadvantage for statistical modeling purposes.

Also, from Equation (2), it can seen that

f(x,−μ,σ)=1σx1−^x2ϕarcsechx−μσ+ϕarcsechx+μσ=f(x,μ,σ).

This means that the pdf shapes of theASHN(μ,σ)distribution coincide with those of theASHN(−μ,σ)distribution. Another remark is that theASHNdistribution can have one mode into(0,1), and it corresponds to the x satisfying the following equation:

2^{σ2 x2}−^σ2arctanh(x)+1−^x2arcsechx=0.

This equation is complex and needs a numerical treatment to determine the value of the mode, if it exists.

The hrf of theASHN(μ,σ)distribution is given by

h(x,μ,σ)=ϕarcsechx+μσ+ϕarcsechx−μσσx1−^x2Φarcsechx+μσ+Φarcsechx−μσ−1.

Some plots off(x,μ,σ)andh(x,μ,σ) are shown in Figure 1.

From Figure 1, the flexibility of the obtained curves is flagrant; J, reversed J, U and bell shapes are observed for the pdf, whereas U, N and reversed J shapes are observed for the hrf. This panel of shapes is a plus for theASHNdistribution, motivating its use for statistical modeling.

3. Distributional Properties

This section is devoted to some mathematical properties satisfied by theASHNdistribution.

3.1. A Likelihood Ratio Order Result

The proposition below shows that theASHNdistribution satisfies a strong intrinsic stochastic order result.

Proposition 2.

LetX∼ASHN(μ,_σ1)andY∼ASHN(μ,_σ2)withμ=0and_σ1>_σ2. Then X is smaller than Y in likelihood ratio order.

In the general case whereμ≠0 , there is no actual proof of such stochastic ordering properties. Further, let us mention that the likelihood order is a strong property, implying various stochastic orders such that the usual stochastic, hazard rate, reversed mean inactivity time, mean residual life and harmonic mean residual life orders, among others. We may refer the reader to [21] for all the theory and details about the concept of stochastic ordering.

3.2. Quantile Function

The theoretical definition of the quantile function (qf) of theASHN(μ,σ) distribution is the inverse function of Equation (1), that is

Q(y,μ,σ)=^F−1(y,μ,σ),y∈(0,1).

In full generality, due to the complexity ofF(x,μ,σ), it is not possible to have a closed-form expression of this qf. However, in the caseμ=0, we arrive at

Q(y,μ,σ)=sechσ^Φ−11−y2,y∈(0,1),

where^Φ−1(x)denotes the inverse function ofΦ(x), which also corresponds to the qf of theN(0,1)distribution. In this case, the first quartile is obtained as_Q1=Q(1/4,μ,σ)≈sechσ×1.150349, the median is given byM=Q(1/2,μ,σ)≈sechσ×0.6744898, and the third quartile is defined by_Q3=Q(3/4,μ,σ)≈sechσ×0.3186394. Further, fromQ(y,μ,σ), one can generate values from theASHNdistribution through basic simulation methods.

3.3. Moments

LetX∼ASHN(μ,σ). As prime definition, for any integer r, by denoting as E the expectation operator, the rth ordinary moment of X is defined by

_mr=E(^Xr)=_∫01 ^xrf(x,μ,σ)dx=1σ_∫01^xr−11−^x2ϕarcsechx+μσ+ϕarcsechx−μσdx.

For the special caseμ=0, one can express it via the qf as

_mr=_∫01 ^{[Q(y,μ,σ)]r}dy=_∫01 ^{sechσΦ−11−y2r}dy.

Clearly, there is no simple expression for_mr. When the parameters are fixed, it can be calculated numerically through standard numerical integration techniques. As the main analytical approach, one can consider a series expansion for_mras stated in the result below.

Proposition 3.

The rth moment ofX∼ASHN(μ,σ)has the following expansion:

_mr=^2r∑k=0+∞−rk^e(2k+r)μM−σ(2k+r),μσ+∑k=0+∞−rk^e−(2k+r)μM−σ(2k+r),−μσ,

whereM(x,a)=E[^exUI(U>a)]withU∼N(0,1),x∈Randa∈R, andI(.)denotes the indicator function.

The functionM(x,a)introduced in Proposition 3 can be viewed as the upper incomplete version of the moment generating function of theN(0,1)distribution. Naturally, it can be bounded from above asM(x,a)≤E(^exU)=^e−x22forx∈Randa∈R. By applying the Markov inequality, a lower is obtained asM(x,a)≥^exa(1−Φ(a))forx≥0anda∈R.

Proposition 3 gives an analytical approach for mathematical manipulations or computations of_mr. Further, the following finite sum approximation is an immediate consequence:

_mr≈^2r∑k=0K−rk^e(2k+r)μM−σ(2k+r),μσ+∑k=0K−rk^e−(2k+r)μM−σ(2k+r),−μσ,

where K denotes a reasonably large integer.

From the moments, we can derive other measures of interest for X. For instance, the mean of X is just_m1, the variance of X can be determined through the Koenig-Huyghens formula involving_m1and_m2, that isV=_m2−_m12, the rth central moment defined by_mrc=E[^(X−_m1)r]can be expressed via_m1,…,_mrby using the binomial formula, the skewness coefficient of X is defined byS=_m3c ^V−32and the kurtosis coefficient of X is given byK=_m4c ^V−2. These coefficients evaluate the “peakedness” and “tailedness” of theASHN distribution, respectively. Figure 2 represents these coefficients while varying the values forμandσ.

From Figure 2, we see that the skewness coefficient can be negative and positive, and the kurtosis coefficient can be either very small or very large. Both have a complex non-monotonic structure. These facts attest to the ability of theASHNdistribution to adapt to various situations from heterogeneous unit data.

3.4. Order Statistics

The order statistics are important since they are involved in many statistical modeling and methods. Here, the basics of them in the context of theASHNdistribution are described. Let_X1,_X2,…,_Xnbe a random sample fromX∼ASHN(μ,σ), and_X(1),_X(2),…,_X(n)be the corresponding order statistics, that is_X(1)≤_X(2)≤…≤_X(n). Then, the pdf of_X(i)has the following general expression:

_fX(i) (x,μ,σ)=n!(i−1)!(n−i)!f(x,μ,σ)^{[F(x,μ,σ)]i−1 [1−F(x,μ,σ)]n−i}.

Owing to Equations (1) and (2), forx∈(0,1), we obtain

_fX(i) (x,μ,σ)=n!(i−1)!(n−i)!1σx1−^x2ϕarcsechx+μσ+ϕarcsechx−μσ×^{2−Φarcsechx+μσ−Φarcsechx−μσi−1}×^{Φarcsechx+μσ+Φarcsechx−μσ−1n−i}.

The pdf of the extreme statistics_X(1)and_X(n)are derived by substitutingi=1andi=nin the above equation, respectively. Other important results are that

E[F(_X(i),μ,σ)]=in+1,V[F(_X(i),μ,σ)]=i(n−i+1)(n+2)⁽ⁿ⁺¹⁾².

The order statistics, as well as their mean and variance, will be useful in the next section. 4. Different Methods of the Parameter Estimation

In this section, we point out some different estimators to estimate the parameters of theASHNmodel. More precisely, the maximum likelihood, maximum product spacings, least squares, weighted least squares, Anderson-Darling and Cramér-von Mises estimates are derived.

4.1. Maximum Likelihood Estimation

Let_X1,_X2,⋯,_Xnbe a random sample from theASHNdistribution with observed values_x1,_x2,⋯,_xn,andΘ=^(μ,σ)Tbe the vector of the model parameters. Then, the log-likelihood function is given by

ℓ=ℓΘ=−nlogσ−n2log2π−∑i=1nlog_xi1−_xi2−12^σ2∑i=1n^{arcsech(_xi)−μ2}+∑i=1nlog1+^{e−2μσ2arcsech(_xi)}.

Based onℓΘ, the maximum likelihood estimations (MLEs) ofμandσ, sayμ^andσ^, respectively, are obtained as

(μ^,σ^)=_{argmaxΘ∈R×(0,+∞)}ℓ(Θ).

Mathematically, this is equivalent to solve the following equations with respect to the parameters:

∂ℓ(Θ)∂μ=1^σ2∑i=1narcsech(_xi)−μ−2^σ2∑i=1narcsech(_xi)^{e−2μσ2arcsech(_xi)}1+^{e−2μσ2arcsech(_xi)}=0

and

∂ℓ(Θ)∂σ=−nσ+1^σ3∑i=1n^{arcsech(_xi)−μ2}+4μ^σ3∑i=1narcsech(_xi)^{e−2μσ2arcsech(_xi)}1+^{e−2μσ2arcsech(_xi)}=0.

From Equation (10), we have

12∑i=1narcsech(_xi)−μ=∑i=1narcsech(_xi)^{e−2μσ2arcsech(_xi)}1+^{e−2μσ2arcsech(_xi)}.

Substituting the right hand side of Equation (12) in Equation (11), the following equation is obtained for the desired solution for^σ2:

^σ2=1n∑i=1n^{arcsech(_xi)−μ2}+2μ∑i=1narcsech(_xi)−μ=1n∑i=1n^{arcsech(_xi)2}−^μ2.

Then, substituting Equation (13) in Equation (9), we obtain the profile log-likelihood according toμas

ℓ(μ)=−n2log1n∑i=1n^{arcsech(_xi)2}−^μ2−n2log2π−∑i=1nlog_xi1−_xi2+∑i=1nlog1+exp−2μarcsech(_xi)^{1n∑i=1narcsech(_xi)2−μ2−1}−∑i=1n^{arcsech(_xi)−μ2}21n∑i=1n^{arcsech(_xi)2}−^μ2.

Following the normal routine of the parameter estimation based on the profile log-likelihood function, we have

∂ℓ(μ)∂μ=nμ1n∑i=1n^{arcsech(_xi)2}−^μ2−μ∑i=1n^{arcsech(_xi)−μ2}−∑i=1narcsech(_xi)−μ1n∑i=1n^{arcsech(_xi)2}−^{μ21n∑i=1narcsech(_xi)2−μ22}−∑i=1n2arcsech(_xi)1n∑i=1n^{arcsech(_xi)2}+4^μ2exp−2μarcsech(_xi)^{1n∑i=1narcsech(_xi)2−μ2−11n∑i=1narcsech(_xi)2−μ22}1+exp−2μarcsech(_xi)^{1n∑i=1narcsech(_xi)2−μ2−1}.

Hence, the numerical methods are needed to obtainμ^. Onceμ^is obtained, the MLEσ^is obtained by taking the square root of^{σ^2} as governed by Equation (13).

The well-known theory of the maximum likelihood method states that, under mild regularity conditions, one can use the bivariate normal distribution with meanμ=(μ,σ)and covariance matrix^I−1, where

I=−_{^∂2∂^μ2ℓ(Θ)^∂2∂μ∂σℓ(Θ)^∂2∂μ∂σℓ(Θ)^∂2∂^σ2ℓ(Θ)Θ=Θ^},

to construct confidence intervals or likelihood ratio test on the parameters. The components of I can be derived through standard derivatives formula. Then, approximate100(1−ϑ)%confidence intervals forμandσcan be determined byμ^±_{zϑ/2 sμ^}andσ^±_{zϑ/2 sσ^}, respectively, where_zϑ/2is the upper(ϑ/2)th percentile of the standard normal distribution,_sμ^is the first diagonal element of^I−1and_sσ^is its second diagonal element. Thus defined, they are the (asymptotic) standard errors (SEs) ofμ^andσ^, respectively.

4.2. Maximum Product Spacing Estimation

Cheng and Aming [22] have proposed the maximum product spacing (MPS) method as an alternative to the maximum likelihood method. It is based on the idea that differences (spacings) between the values of the cdf at consecutive data points should be identically distributed. Now, let_X(1),_X(2),…,_X(n)be the order statistics from theASHNdistribution with sample size n, and_x(1),_x(2),…,_x(n)be the ordered observed values. Then, the MPS estimates (MPSEs) ofμandσ, say_μ^MPSand_σ^MPS, respectively, are given as

(_μ^MPS,_σ^MPS)=_{argmaxΘ∈R×(0,+∞)}MPS(Θ),

where

MPS(Θ)=1n+1∑i=1n+1logF(_x(i),μ,σ)−F(_x(i−1),μ,σ).

They are also given as the simultaneous solutions of the following equations:

∂MPS(Θ)∂μ=1n+1∑i=1n+1_Fμ^′ (_x(i),μ,σ)−_Fμ^′ (_x(i−1),μ,σ)F(_x(i),μ,σ)−F(_x(i−1),μ,σ)=0

and

∂MPS(Θ)∂σ=1n+1∑i=1n+1_Fσ^′ (_x(i),μ,σ)−_Fσ^′ (_x(i−1),μ,σ)F(_x(i),μ,σ)−F(_x(i−1),μ,σ)=0,

where

_Fμ^′ (x,μ,σ)=1σϕarcsechx−μσ−ϕarcsechx+μσ

and

_Fσ^′ (x,μ,σ)=1^σ2arcsechx−μϕarcsechx−μσ+arcsechx+μϕarcsechx+μσ.

4.3. Least Squares Estimation

The least square estimates (LSEs) ofμandσ, say_μ^LSEand_σ^LSE, respectively, are obtained as

(_μ^LSE,_σ^LSE)=_{argminΘ∈R×(0,+∞)}LSE(Θ),

where

LSE(Θ)=∑i=1n^{F(_x(i),μ,σ)−EF(_X(i),μ,σ)2},

where, by Equation (8),EF(_X(i),μ,σ)=i/(n+1)fori=1,2,…,n. Then,_μ^LSEand_σ^LSEare solutions of the following equations:

∂LSE(Θ)∂μ=2∑i=1n_Fμ^′ (_x(i),μ,σ)F(_x(i),μ,σ)−in+1=0

and

∂LSE(Θ)∂σ=2∑i=1n_Fσ^′ (_x(i),μ,σ)F(_x(i),μ,σ)−in+1=0,

where_Fμ^′ (_x(i),μ,σ)and_Fσ^′ (_x(i),μ,σ)are mentioned before.

4.4. Weighted Least Squares Estimation

Similarly to LSEs, the weighted least square estimates (WLSEs) ofμandσ, say_μ^WLSEand_σ^WLSE, respectively, are given as

(_μ^WLSE,_σ^WLSE)=_{argminΘ∈R×(0,+∞)}WLSE(Θ),

where

WLSE(Θ)=∑i=1n1VF(_X(i),μ,σ)^{F(_x(i),μ,σ)−EF(_X(i),μ,σ)2},

where, by Equation (8),EF(_X(i),μ,σ)=i/(n+1)andVF(_X(i),μ,σ)=i(n−i+1)/[(n+2)⁽ⁿ⁺¹⁾²]fori=1,2,…,n. Then,_μ^WLSEand_σ^WLSEare solutions of the following equations:

∂WLSE(Θ)∂μ=2∑i=1n(n+2)⁽ⁿ⁺¹⁾²i(n−i+1)F(_x(i),μ,σ)−in+1_Fμ^′ (_x(i),μ,σ)=0

and

∂WLSE(Θ)∂σ=2∑i=1n(n+2)⁽ⁿ⁺¹⁾²i(n−i+1)F(_x(i),μ,σ)−in+1_Fσ^′ (_x(i),μ,σ)=0.

4.5. Anderson-Darling Estimation

The Anderson-Darling minimum distance estimates (ADEs) ofμandσ, say_μ^ADand_σ^AD, respectively, are determined as

(_μ^AD,_σ^AD)=_{argminΘ∈R×(0,+∞)}AD(Θ),

where

ADΘ=−n−∑i=1n2i−1nlog[F(_x(i),μ,σ)]+log1−F(_x(n+1−i),μ,σ).

Therefore,_μ^ADand_σ^ADcan be obtained as the solutions of the following system of equations:

∂ADΘ∂μ=−∑i=1n2i−1n_Fμ′(_x(i),μ,σ)F(_x(i),μ,σ)−_Fμ′(_x(n+1−i),μ,σ)1−F(_x(n+1−i),μ,σ)=0

and

∂ADΘ∂σ=−∑i=1n2i−1n_Fσ′(_x(i),μ,σ)F(_x(i),μ,σ)−_Fσ′(_x(n+1−i),μ,σ)1−F(_x(n+1−i),μ,σ)=0.

4.6. The Cramér-von Mises Estimation

The Cramér-von Mises minimum distance estimates (CVMEs) ofμandσ, say_μ^CVMand_σ^CVM, respectively, are specified as

(_μ^CVM,_σ^CVM)=_{argminΘ∈R×(0,+∞)}CVM(Θ),

where

CVMΘ=112n+∑i=1n^{F(_x(i),μ,σ)−2i−12n2}.

Therefore, the estimates_μ^CVMand_σ^CVMcan be obtained as the solutions of the following system of equations:

∂CVMΘ∂μ=2∑i=1nF(_x(i),μ,σ)−2i−12n_Fμ′(_x(i),μ,σ)=0

and

∂CVMΘ∂σ=2∑i=1nF(_x(i),μ,σ)−2i−12n_Fσ′(_x(i),μ,σ)=0.

All the presented equations contain complex non-linear functions; it is not possible to obtain explicit forms of all estimates. Therefore, they need to be solved through numerical methods such as the Newton-Raphson and quasi-Newton algorithms. In addition, Equations (9) and (14)–(18) can be also optimized directly by using the software such as R (constrOptim and optim), S-Plus and Matlab to numerically optimizeℓ(Θ),MPSΘ,LSEΘ,WLSEΘ,ADΘandCVMΘfunctions.

5. Empirical Simulations

In this section, we perform two graphical simulation studies to see the performance of the above estimates with varying sample size n. We generateN=1000samples of sizen=20,25,…,1000from theASHNdistribution based on the following parameter values: (μ=2,σ=2) and (μ=0.5,σ=0.5 ) for the first and second simulation studies, respectively. The random numbers generation is obtained by the qf of the model. All the estimates based on the estimation methods are obtained by using the constrOptim function in the R program. Further, we calculate the empirical mean, bias and mean square error (MSE) of the estimates for comparisons between the methods. Forϵ=μorϵ=σ, the bias and MSE associated toϵare calculated by

Bia_sϵ(n)=1N∑i=1N(ϵ−_ϵ^i),MS_Eϵ(n)=1N∑i=1N^{(ϵ−_ϵ^i)2},

respectively, where i is related to the ith sample. We expect that the empirical means are close to true values when the MSEs and biases are near zero. The results of this simulation study are shown in Figure 3 and Figure 4.

Figure 3 and Figure 4 show that all estimates are consistent since the MSE and biasedness decrease to zero with increasing sample size as expected. One can state that all estimates are asymptotic unbiased. According to these two simulation studies, the amount of the biases and MSEs of the MLE method are smaller than those of the other methods for both parameters. Therefore, the ;MLE method can be chosen as more reliable than other methods of the newly defined model. Generally, the performances of all estimates are close when sample size increases. The similar results can be seen for different parameter values.

Moreover, we also give simulation study of the MLEs based on their 95% confidence intervals. In this regard, we use the coverage probability (CP) criteria defined by

C_Pϵ(n)=1N∑i=1NI(_ϵi^±1.95996_sϵi^),

where_sϵ^i is the SE of the MLE_ϵ^i . Figure 5 displays the obtained simulation results. From Figure 5, as expected, for each parameter, the CPs converge to the nominal value, that is 0.95, when sample size increases. The simulation results verify the consistency property of the MLEs.

6. A New Quantile Regression Model Based on the SpecialASHNDistribution 6.1. Motivation

The quantile regression has been developed in the seminal work of [19] as a way to model the conditional quantiles of an outcome variable as a function of covariates (regressors). Since this analysis aims to model the conditional quantiles of the response variable, it is a good robust alternative model to the ordinary LSE model, which estimates the conditional mean of the response variable. This is because the mean is affected by a skewed distribution or outliers in the measurements. Hence, the quantile response regression model will be less sensitive to outliers than the mean response regression model.

On the other hand, if the support of the response variable is defined on the unit interval, one can use an unit regression model based on an unit distribution for modeling the conditional mean or quantiles of the response variable via covariates. The beta regression [23] model first comes to mind to relate to continuous unit mean response variables in the unit interval with covariates. One may also see [12,13,24,25,26] for alternative unit mean response regression models to beta regression models. If the conditional dependent variable is skewed or has outliers, the quantile response modeling may be more appropriate when compared with the mean response modeling. The model is also motivated by the natural idea of replacing mean by median as a central tendency measure when the response data is severely asymmetric [27].

On the other hand, with the re-parameterizing the probability distribution as a function of the quantile approach, the Kumaraswamy [28,29] and unit Weibull [30] quantile regression models have been proposed for modeling the conditional quantiles of the unit response. One may also refer to [27,31,32,33,34] for alternative quantile response regression models. On the basis of these references, we want to propose an alternative quantile regression model considering a parameterization of theASHNdistribution in terms of its any quantile. More precisely, the re-parameterizing process is applied via a scale parameter as being a quantile of the of theASHNdistribution.

6.2. Proposed Quantile Regression Model

Now, we can focus on introducing an alternative quantile regression model based on a specialASHN(μ,σ)distribution. Since theASHNdistribution has not an explicit qf, we propose another distribution based on a specialASHNdistribution. We call it exponentiatedASHN(EASHN) distribution. Its cdf and pdf are given by

G(y,α,σ)=^{2−2Φ1σarcsechyα}

and

g(y,α,σ)=2ασy1−^y2ϕ1σarcsechy^{2−2Φ1σarcsechyα−1},

respectively, wherey∈(0,1)andα,σ>0. Fory∉(0,1) , standard completions on these functions are performed. The cdf in Equation (19) is obtained as the exponentiatedASHN(0,σ)distribution, that isG(y,α,σ)=^F(y,0,σ)α. We can call this model as Lehmann type IASHNmodel.

The qf of theEASHNdistribution is given by

Q(u,α,σ)=sechσ^Φ−11−^u1α2,

whereu∈(0,1). We are also motivated with the quantile regression modeling thanks to its manageable qf. Then, the pdf of theEASHNdistribution can be re-parameterized in terms of its uth quantile asη=Q(u,α,σ). Letσ=arcsechη/^Φ−1(2−^u1/α)/2. Then, the cdf and pdf of the re-parameterized distribution are given by

G(y,α,η)=^{2−2ΦΦ−11−u1α/2arcsechηarcsechyα}

and

g(y,α,η)=2α^Φ−11−^u1α/2arcsechηy1−^y2ϕ^Φ−11−^u1α/2arcsechηarcsechy×^{2−2ΦΦ−11−u1α/2arcsechηarcsechyα−1},

respectively, whereα>0is the shape parameter. The parameterη∈(0,1) represents the quantile parameter and it is assumed that u is known. A random variable Y having the pdf in Equation (21) is denoted byY∼EASHN(α,η,u) . Some possible shapes of the re-parameterized model are shown in Figure 6. We see that the possible pdf shapes of theEASHNdistribution are the skewed shapes as well as U-shapes, N-shapes and increasing shapes.

We present the quantile regression model based on theEASHN distribution with pdf in Equation (21). Let_y1,_y2,…,_ynbe n random observations from the re-parameterized distribution such that, fori=1,…,n,_yiis a realization ofY∼EASHN(α,_ηi,u), with unknown parameters_ηiandβ, recalling that the parameter u is known. Then theEASHNquantile regression model is defined as

g(_ηi)=_xi ^βT,

whereβ=^{_β0,_β1,_β2,…,_βpT}and_xi=1,_xi1,_xi2,_xi3,…,_xipare the unknown regression parameter vector and known ith vector of the covariates. Thus defined,g(x)is the link function which is used to link the covariates to conditional quantile of the response variable. For instance, whenu=0.5, the covariates are linked to conditional median of the response variable. The choice of the appropriate link function should be done considering the domain of the distribution.

6.3. Parameter Estimation

The unknown parameters of theEASHNquantile regression model are obtained by means of the MLE method. Since theEASHNdistribution is defined on the unit interval, we use the logit-link function, that is

g(_ηi)=logit(_ηi)=log_ηi1−_ηi=_xi ^βT,

implying that

_ηi=exp_xi ^βT1+exp_xi ^βT.

By putting Equation (22) into Equation (21), the log-likelihood function of theEASHNquantile regression model is

ℓΩ=nlog2−n2log2π+nlogα+nlog^Φ−11−^u1α2−∑i=1nlogarcsech(_ηi)_yi1−_yi2−12^{Φ−11−u1α22}∑i=1n^{arcsech(_yi)arcsech(_ηi)2}+(α−1)∑i=1nlog2−2Φarcsech(_yi)arcsech(_ηi)^Φ−11−^u1α2,

whereΩ=α,β denotes the unknown parameter vector. Since Equation (23) includes nonlinear function according to model parameters, it can be maximized directly by software such as R, S-Plus, and Mathematica. Note that, whenu=0.5, this is equivalent to modeling the conditional median. Under mild conditions of regularity, the asymptotic distribution ofΩ^−Ωis multivariate normal_Np+10,^J−1, where variance-covariance matrix^J−1 is defined by the inverse of the expected information matrix. For practical aims, we can use the observed information matrix instead of J. The elements of this observed information matrix are evaluated numerically by the software. We use the maxLik function implemented in the R software to maximize Equation (23) (see [35]). This function also gives asymptotic SEs numerically, which are obtained by the observed information matrix.

6.4. Residual Analysis

Residual analysis may be necessary to verify if the regression model is suitable. To see this, we work with the randomized quantile residuals [36] and the Cox-Snell residuals [37].

Fori=1,…,n, the ith randomized quantile residual is defined by

_r^i=^Φ−1G(_yi,α^,_η^i),

whereG(y,α,η)is the cdf of the re-parameterizedEASHN distribution specified by Equation (20) and_η^i is defined by Equation (22) withβ^replacingβ. If the fitted model successfully processes the dataset, the distribution of the randomized quantile residuals will distribute the standard normal distribution.

Alternatively, fori=1,…,n, the ith Cox and Snell residual is given by

_e^i=−log1−G(_yi,α^,_η^i).

If the model fits to data accordingly, the distribution of the Cox and Snell residuals will distribute a standard exponential distribution, that is with scale parameter 1. 7. Data Analysis

To emphasize the importance of the modeling ability of theASHNnormal distribution, this section is devoted to three real data applications for both univariate data and quantile modelings.

7.1. Univariate Real Data Modeling

Here, we provide applications to two real datasets to prove empirically the potentiality of theASHNmodel. The proposed model is compared with some well-known two-parameter unit distributions in the literature, namely:

• Beta distribution.

  The two-parameter beta pdf is given by

  _fBeta(x,μ,σ)=1Bμ,σ^{xμ−1 1−xσ−1},x∈(0,1),

  whereμ>0andσ>0are shape parameters, andB(μ,σ)is the classical beta function.

• Kumaraswamy (Kw) distribution (see [3]).

  The two-parameter Kw pdf is expressed as

  _fKw(x,μ,σ)=μσ^{xμ−1 1−xμσ−1},x∈(0,1),

  whereμ>0andσ>0are shape parameters.

• Johnson_SB distribution (see [1]).

  The two-parameter Johnson_SBpdf is given by

  _fSB (x,μ,σ)=σx1−xϕσlogx1−x+μ,x∈(0,1),

  whereμ∈Randσ>0are shape parameters, andϕ(·)is the pdf of the standard normal distribution. For each model, we estimate the unknown parameters using the maximum likelihood approach.

Two datasets are considered. For them, in order to determine the optimum model, we compute the estimated log-likelihood valuesℓ^, Akaike Information Criteria (AIC), Bayesian Information Criteria (BIC), Kolmogorov-Smirnov (KS), Anderson-Darling (^A*) and Cramér-von Mises (^W*) goodness-of-fit statistics for all models. In general, it can be chosen as the best model the one with the smaller values of the AIC, BIC, KS,^A*and^W*statistics and the larger values ofℓ^ . The p-value of the KS test is also considered; more it is close to 1, better is the model. All computations are performed by using the maxLik [35] and goftest routines in the R software.

7.1.1. Data Analysis I

First, we consider an application to real dataset to show the modeling ability of the proposed distribution. The dataset introduces failure times of the 20 mechanical components given in [38]. The data are: 0.067, 0.068, 0.076, 0.081, 0.084, 0.085, 0.085, 0.086, 0.089, 0.098, 0.098, 0.114, 0.114, 0.115, 0.121, 0.125, 0.131, 0.149, 0.160, 0.485. Recently, these data was analyzed via different approaches by [15,39,40].

Table 1 lists the MLEs of the parameters and their SEs from the above fitted models and the values of the statistics:ℓ^, AIC, BIC,^A*,^W*and KS goodness-of-fit statistics. As it can be seen, the results indicate that theASHNmodel has the smallest values of these statistics among the fitted models, and therefore it could be considered as the best model. The p-value of the KS test confirms this claim.

The plots of the fitted pdfs and cdfs are displayed in Figure 7. These plots show that theASHNmodel provides the correct fit to these data compared to other models. Further, it captures data skewness and kurtosis better than other models.

Figure 8 shows plots of the profile log-likelihood (PLL) functions for the parametersμandσbased on the first dataset. We observe that the likelihood equations have unique solutions for the MLEs.

7.1.2. Data Analysis II

Here, the Better Life Index (BLI) dataset is used to illustrate the usefulness of theASHN distribution. The dataset can be found via link https://stats.oecd.org/index.aspx?DataSetCode=BLI2015. The BLI dataset is used to classify the OECD (Organisation for Economic Co-operation and Development) countries with 11 indicators and 24 variables as well as non-OECD economies such as Brazil and Russia. Here, we use an indicator that is entitled Job security as the dataset. This indicator presents the probability to become unemployed. Recently, these data was analyzed by [14]. We give the summary statistics of the dataset in Table 2. The data are right-skewed and have a consequent kurtosis.

Table 3 lists the MLEs, their SEs,ℓ^ and goodness-of-fits statistics from the fitted models for this dataset. Table 3 shows that the proposed model could be chosen as the best model among the fitted models since it has the lowest values of the AIC, BIC,^A*,^W*and KS statistics and have the biggestℓ^value. It also has the biggest p-value of the KS test.

The plots of the fitted pdfs and cdfs are displayed in Figure 9. These plots show that theASHNmodel provides the good fit to these data compared to the other models.

Figure 10 shows plots of PLL functions for the parametersμandσbased on the second dataset. From this figure, we see that the likelihood equations have unique solutions for the MLEs.

7.2. The Quantile Modeling Application of the Reading Accuracy with the Dyslexia and Intelligence Quotient

Here, a real data application is given in order to see the applicability of the newly defined regression model. We compare its results with the unit Weibull quantile regression model [30]. The quantile parameter u has been taken as 0.5 to model the median for regression models. The pdf of the unit Weibull quantile distribution is given by

f(y,α,μ)=αlog0.5ylogη^{logylogηα−1}0.^{5logy/logηα} ,y∈(0,1),

whereη∈(0,1)is the median andα>0 is the shape parameter. The dataset consists of the reading accuracy for nondyslexic and dyslexic Australian children and contains 44 observations on 3 variables. The variable of interest is accuracy providing the scores on a test of reading accuracy taken by 44 children, which is predicted by the two regressors: dyslexia and nonverbal intelligence quotient (IQ). The dataset has been collected by [41], and analyzed by [42,43] via the beta regression modeling based on the data mean modeling. It is noticed that the original reading accuracy score has been transformed by [43] so that accuracy is in the open unit interval. Further, this dataset can be found easily via betareg function [42] in the R software.

The aim is to associate the reading accuracy values(y)with covariates. The response variable and covariates are:

• y: reading score;

• _x1: Is the child dyslexic? (0 for no, 1 for yes);

• _x2: nonverbal intelligence quotient (IQ, converted to z scores).

The regression model based on_ηiis given by

logit(_ηi)=_β0+_{β1 xi1}+_{β2 xi2},i=1,2,…,44,

where_ηidenotes the median for the unit Weibull andEASHNmodels.

The results of theEASHN and unit Weibull regression models with model selection criteria are given in Table 4. As seen from the values of AIC and BIC statistics, the proposed regression model has lower values than those of the unit Weibull regression model. So, one can say that theEASHNregression model exhibits better modeling ability than the unit Weibull regression model. Additionally, according to the estimated parameters of theEASHNregression model, the parameters_β1and_β2have been seen statistically significant at any usual level. Hence, it is concluded that, when IQ increases, the reading accuracy increases also. However, the reading accuracy of the children with no dyslexia is higher than those of the children with dyslexia as expected.

Figure 11 and Figure 12 display the QQ plots of the randomized quantile residuals and PP plot of the Cox-Snell residuals for both regression models, respectively. These figures indicate that the fit of theEASHNregression model is better than the one of the unit Weibull model.

Since the randomized quantile residuals have standard normal distributions, one may see whether they fit this corresponding distribution. The KS,^A*and^W* results are given in Table 5. It is clear that the results based on theEASHNquantile regression model of the randomized quantile residuals are more suitable than those of the unit Weibull regression model.

8. Conclusions We define a new unit model, called “arcsech” normal distribution, in order to model percentage, proportion and rate measurements. The idea is to take advantage of the hyperbolic arcsecant function to transpose the modeling capacities of the normal distribution for the processing of data defined on the unit interval. We investigate general structural properties of the new distribution. The model parameters are estimated by six different methods. The simulation studies are performed to see the performances of these estimates. The empirical findings indicate that the proposed model provides better fits than the well-known unit probability distributions in the literature for both its univariate data modeling and its regression modeling. It is hoped that the new distribution will attract attention in the other disciplines.

Model	μ^	σ^	ℓ^	AIC	BIC	A*	W*	KS
ASHN	2.9179	0.4322	33.2443	−62.4885	−60.4970	1.1850	0.1664	0.1746
	(0.0966)	(0.0684)						[0.5754]
Beta	3.1126	21.8245	27.8813	−51.7626	−49.7711	2.2611	0.3726	0.2537
	(1.0287)	(7.7997)						[0.1521]
Kw	1.5877	21.8673	25.6484	−47.2968	−45.3054	2.6889	0.4681	0.2626
	(0.3966)	17.9755						[0.1265]
JohnsonSB	3.8952	1.8605	31.3599	−58.7198	−56.7283	1.5531	0.2307	0.2039
	(0.6554)	(0.2942)						[0.3765]

Minimum	Mean	Median	Maximum	Variance	Skewness	Kurtosis	n
0.0240	0.0567	0.0515	0.1780	0.0007	2.7117	12.0173	36

Model	μ^	σ^	ℓ^	AIC	BIC	A*	W*	KS
ASHN	3.6422	0.3791	90.1076	−176.2152	−173.0481	0.5963	0.0895	0.1261
	(0.0632)	(0.0447)						[0.6162]
Beta	5.8569	97.1458	86.9760	−169.9519	−166.7848	1.1152	0.1768	0.1636
	(0.5166)	(6.2564)						[0.2903]
Kw	2.1577	373.3878	82.0487	−160.0975	−156.9305	2.2041	0.3651	0.1916
	(0.0648)	8.4525						[0.1422]
JohnsonSB	7.1149	2.4608	89.6573	−175.3146	−172.1476	0.6666	0.1008	0.1322
	(0.8440)	(0.2864)						[0.5554]

Parameters	EASHN			Unit-Weibull
	Estimate	SE	p-Value	Estimate	SE	p-Value
β0	2.2810	0.0025	<0.001	2.4045	0.2589	<0.001
β1	−1.0490	0.0028	<0.001	−1.3362	0.3751	0.0003
β2	0.5918	0.00001	<0.001	0.4837	0.2453	0.0486
α	0.1260	0.00001	<0.001	0.9795	0.1193	<0.001
ℓ	37.9466			37.3185
AIC	−67.8934			−66.6369
BIC	−60.7566			−59.5001

Models	KS	p-Value	A*	p-Value	W*	p-Value
EASHN	0.0849	0.9093	0.4211	0.8267	0.0502	0.8775
Unit-Weibull	0.1159	0.5955	0.4989	0.7470	0.0720	0.7419

Author Contributions

M.Ç.K., C.C. and Z.S.K. have contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

The authors are grateful to the two anonymous referees helpful comments that improved the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The proofs of our main results are contained in this appendix.

Proof of Proposition 1.

Firstly, it is noticed that the hyperbolic secant function is the symmetrical on the-∞,∞interval as well as it is the increasing function on the-∞,0interval and it is the decreasing function on the0,∞interval. Based on the representationX=sechY, whereY∼N(μ,^σ2), the cdf of X can be determined as

F(x)=PX≤x=P-∞≤Y≤-arcsechx+Parcsechx≤Y≤∞=Φ-arcsechx-μσ+1-Φarcsechx-μσ=2-Φarcsechx+μσ-Φarcsechx-μσ.

We get the declared definition ofF(x,μ,σ). By differentiation ofF(x,μ,σ)with respect to x, since∂(arcsechx)/∂x=-^x1-x2-1, the pdff(x,μ,σ)follows, ending the proof. □

Proof of Proposition 2.

In the caseμ=0 , owing to Equation (4), some simplifications give

f(x,0,_σ1)f(x,0,_σ2)=2ϕ(arcsechx)/_σ1/(_σ1x1-^x2)2ϕ(arcsechx)/_σ2/(_σ2x1-^x2)=_{σ2 σ1}^{e12(_σ1-_σ2)_σ1+_{σ2σ12 σ22}(arcsechx)2}.

Sincearcsechxis a positive decreasing function, if_σ1>_σ2, the above ratio function is decreasing with respect to x as composition of an increasing exponential function and a decreasing function. This proves the desired result. □

Proof of Proposition 3.

We propose to exploit the following characterization of theASHNdistribution: We can express X asX=sechYwithY∼N(μ,^σ2). Then, through the use of the generalized version of the binomial formula, we obtain

^Xr=^(sechY)r=^{2rerY (1+e2Y)r}I(Y<0)+I(Y=0)+^{2re-rY (1+e-2Y)r}I(Y>0)=^2r∑k=0+∞-rk^e(2k+r)YI(Y<0)+I(Y=0)+^2r∑k=0+∞-rk^e-(2k+r)YI(Y>0).

Therefore, sinceP(Y=0)=0, we get

_mr=^2r∑k=0+∞-rkE[^e(2k+r)YI(Y<0)]+∑k=0+∞-rkE[^e-(2k+r)YI(Y>0)].

In the distribution sense, one can writeY=μ+σUwithU∼N(0,1), implying that

E[^e(2k+r)YI(Y<0)]=^e(2k+r)μE[^eσ(2k+r)UI(U<-μ/σ)]=^e(2k+r)μM-σ(2k+r),μσ

and

E[^e-(2k+r)YI(Y>0)]=^e-(2k+r)μM-σ(2k+r),-μσ.

The stated result follows by combining the equations above together. □