# Fully Bayesian Estimation of Simultaneous Regression Quantiles under Asymmetric Laplace Distribution Specification

Authors: Josephine Merhi Bleik
Journal: Journal of Probability and Statistics
 (2019)

---

## Abstract

In this paper, we are interested in estimating several quantiles simultaneously in a regression context via the Bayesian approach. Assuming that the error term has an asymmetric Laplace distribution and using the relation between two distinct quantiles of this distribution, we propose a simple fully Bayesian method that satisfies the noncrossing property of quantiles. For implementation, we use Metropolis-Hastings within Gibbs algorithm to sample unknown parameters from their full conditional distribution. The performance and the competitiveness of the underlying method with other alternatives are shown in simulated examples.


---

## Body

## 1. Introduction

Quantile regression has received increasing attention from both theoretical and empirical points of view. It is a statistical procedure intended to estimate conditional quantile of a response distribution, given a set of covariates, and to conduct inference about them.In regression, when we deal with highly skewed data distribution, traditional mean regression may not explore interesting aspects of the relationship between the response variable and the available covariates. Furthermore, in the presence of outliers, the evaluation of the response average becomes much more complicated. However, since a set of quantiles often provides a complete description of the response distribution conditionally to the given covariates, quantile regression offers a practical alternative to the traditional mean regression.Quantile regression was initially introduced by Koenker and Bassett [1], where inference proceeds by minimizing a check loss function. Based on their theory, many other approaches have been emerged to study the regression quantiles in both frequentist and Bayesian points of view.In the frequentist framework, different techniques have been developed in the literature to infer regression quantiles. To estimate linear regression quantiles, Koenker and Bassett [1] used linear programming technique, whereas Gutenbrunner and Jurecková [2] considered the corresponding dual problem whose solution represents the quantile coefficients estimator. From theoretical point of view, asymptotic normality and consistency properties are proven in Koenker and Bassett [1], Zhou et al. [3], and Koenker and Xiao [4]. For nonparametric quantile regression, spline methods are commonly used for quantile regression fitting (see Nychka et al. [5], Koenker and Portnoy [6], and Koenker et al. [7]). Oh et al. [8] proposed, in univariate and bivariate settings, an approximation of the nonlinear optimization problem by a sequence of penalized least squares type nonparametric mean estimation problems. Koenker and Mizera [9] explored the estimation of triograms in bivariate setting. Koenker et al. [10] proposed an additive models estimation method for quantile regression based on selecting the smoothing parameter and constructing confidence bands for nonparametric components, in univariate and bivariate settings. The univariate regression quantile spline is extended to problems with several covariates in He and Ng [11] and Li et al. [12] where, in this latter, regression quantiles are considered in reproducing kernel Hilbert spaces (RKHS).In the Bayesian framework, different methods have been successfully emerged to study quantile regression. This was pioneered by Yu and Moyeed [13], where the authors specified the Bayesian quantile regression model and assumed that the error terms are iid according to the asymmetric Laplace distribution (ALD), denoted by ALD(μ,σ=1,p) with μ∈R being its p-th quantile, σ is the scale parameter, and p∈[0,1] are the asymmetry parameters. They also showed that the Bayesian approach yields a proper full conditional distribution, even for an improper uniform prior on the quantile coefficients and they used the Metropolis-Hastings algorithm for implementation. Further, Kozumi and Kobayashi [14] and Bernardi et al. [15] considered that the scale parameter σ of the ALD is unknown. They used the location-scale mixture representation of the ALD, which enables developing a Gibbs sampling algorithm for the Bayesian quantile regression model. In the literature, we can also find other methods that extended quantile regression model. For censored data, Kobayashi et al. [16] proposed Tobit quantile regression with endogenous variables and Li and Bradic [17] developed regression adjustment for local adaptive representation of random forests. In continuous dependent variable case and for variable selection and coefficient estimation, Alhamzawi et al. [18] used adaptive LASSO quantile regression and, recently, Sayed-Ahmed [19] applied this approach for small sample sizes, whereas Abbas and Thaher [20] developed Bayesian adaptive Tobit regression. Benoit and Van den Poel [21] proposed Bayesian quantile regression methods for binary response data and Alhamzawi and Ali [22] adapted the quantile regression model to deal with longitudinal ordinal data.All these methods mentioned above, from both frequentist and Bayesian points of view, studied single regression quantile. However, there are many applications where tackling several quantiles at the same time is needed. For this direction, if one uses methods that deal with single quantile regression to estimate regression quantiles separately, results may not be satisfactory. Indeed, it is known that the separate estimation of a set of quantiles may break the quantile monotonicity property (see He [23]), which means that quantile curves may cross.Fortunately, many approaches have been proposed to overcome the crossing problem. In the frequentist framework, Liu and Wu [24] and Bondell et al. [25] considered monotonicity constraints for a fixed number of quantiles in both linear and nonparametric cases. Sangnier et al. [26] proposed a nonparametric method that estimates multiple conditional quantiles simultaneously when assuming that the quantile functions spans an RKHS. In this work, they choose an appropriate matrix-valued kernel taking into account the distance between quantiles.On the other hand, in simultaneous Bayesian framework, Reich et al. [27] proposed a two-stage semiparametric method for linear quantile regression, which is one of the first methods addressing the crossing problem of quantiles. At the first stage, they used Koenker and Bassett’s [1] method to estimate quantile coefficients. At the second stage, these coefficients are reestimated with Bernstein polynomial basis function with constraints on their coefficients. Later, Reich [28] proposed a fully Bayesian approach based on interpolating the quantile function on different quantile levels using a piecewise Gaussian basis function. Tokdar et al. [29] proposed a simultaneous linear quantile regression method for a univariate explanatory variable with a bounded and convex domain. They showed that a linear quantile, qτ(x)=β0(τ)+xβ(τ), is monotonically increasing if and only if β0(τ) and β(τ) are linear combination of two monotonic increasing functions in τ∈[0,1]. Afterwards, Yang and Tokdar [30] extended this method to the multivariate explanatory case. Another two-stage method is proposed by Rodrigues and Fan [31]. At the first stage, they used the standard Bayesian quantile regression of Yu and Moyeed [13] fitted separately at the different quantile levels. They introduced induced quantiles, at the second stage, using a relation between quantiles of the ALD distribution. At the same stage, they used the Gaussian process regression that monotonizes quantile functions by borrowing strength from the nearby ones. However, this method is still affected by the outputs of the first stage since different likelihoods are used for initial estimates. For nonparametric quantile regression, a cubic B-spline method is considered by Muggeo et al. [32], solving an L2 penalisation problem with monotonicity constraint on spline coefficients. Many recent works on simultaneous Bayesian quantile regression can be found in literature such as Das and Ghosal [33], Das and Ghosal [34], Rodrigues et al. [35], and Petrella and Raponi [36].In this paper, our contribution aims at developing a fully Bayesian approach that estimates multiple quantiles simultaneously in one step. Unlike the other Bayesian methods, our proposed method is not based on misspecification, i.e., the ALD is the real underlying response distribution. We used a relation between quantiles of the ALD distribution and Metropolis-Hastings within Gibbs algorithm for implementation. Our proposed method provides parallel quantile functions estimators and, thus, the noncrossing is a natural result. Our numerical results show, using specific empirical criteria, that the estimated regression quantiles respect the quantile level ordering.The paper is organised as follows. Section2 describes the method presenting the model and establishing the relation between quantiles of the ALD distribution used in the estimation procedure. The algorithm associated with the proposed method is given in Section 3, where we introduce the location-scale mixture representation of the ALD. In Section 4, simulated examples are provided to show the performance and the flexibility of the method, and the final section presents our conclusion.


## 2. Simultaneous Bayesian Estimation of Quantiles

### 2.1. Model

We consider the quantile regression model(1)Yi=qpXi+ϵi,i=1,…,n,with (ϵi)i~iidALD(μ=0,σ,p), where μ, σ, and p, respectively, are the location, the scale, and the asymmetry parameters. This leads to assuming that the response variable Yi, given the covariable Xi, follows an ALD distribution:(2)Yi∣qpXi,σ,p~ALDqpXi,σ,p∀i=1,…,n,whose probability density function (p.d.f.) is given by(3)fyi;qpxi,σ,p=p1-pσexp⁡-ρpyi-qpxiσ,where ρp(u)=u(p-1u<0) is the check loss function. The explanatory variables Xi’s are supposed to be iid distributed according to an arbitrary continuous distribution, PX, whose support is Xd, X⊂R, d≥1. In what follows, the function qp(.) and the parameters σ, p, and PX are supposed to be unknown.In all subsequent sections, any quantity in bold represents a vector; capital letters denote random variables or vectors, whereas lowercase letters denote observed values of the corresponding random vectors.


### 2.2. Bayesian Procedure

#### 2.2.1. Likelihood

To infer simultaneouslys distinct quantiles with a fully Bayesian approach, say qτ1,…,qτs with distinct orders τ1,…,τs, of the conditional distribution of Y∣X, one needs to characterize the likelihood through all these s unknown quantiles.As explained below, this is done first by partitioning the whole sample intos subsamples, which requires a sufficient number of observations, second by using the relation between any two quantiles of the ALD distribution, and third by rewriting the whole likelihood through s terms, where the j-th term only depends on qτj.(1)
Consider the following partition of the whole sample throughs subsamples:(4)XIj,YIj=Xi,Yii∈Ij,∀j∈1,…,s,whereIj=(j-1)r+1,…,jr with Ij⋂Ik=∅, ∀j≠k, and j,k∈{1,…,s}; suppose that r=n/s is integer, without loss of generality.(2)
To characterize the likelihood through all quantiles of interest, we use the relation that links anyτ-th quantile to the p-th quantile of the underlying ALD distribution (see Yu and Zhang [37] and Rodrigues and Fan [31]):(5)qτX=qpX+σgτ,p,whereg(τ,p)=1/1-plog⁡τ/p10<τ≤p-1/plog⁡(1-τ/1-p)1p<τ<1.(3)
From (5), we rewrite the model given by (1) as follows:(6)YIj=qτj,Ij-σgτj,p+ϵIj,∀j=1,…,s,where(7)YIj=Yijij∈Ij,qτj,Ij=qτjXijij∈Ij,gτj,p=gτj,p1ij∈Ijij∈Ij,ϵIj=ϵijij∈Ij.Then, the likelihood associated with the model given by (1) is rewritten as the product of s likelihoods; each one only depends on one subsample: (8)Lqτ1,I1,…,qτs,Is,σ,p;xI1,…,xIs,yI1,…,yIs=∏j=1s∏ij∈Ijp1-pσexp⁡-ρpyij-qτjxij-σgτj,pσ.


#### 2.2.2. Priors

In this section, we specify the prior distributions for all unknown quantities. It is worth mentioning thatPX is unknown even it is not our primary interest, and therefore it should require a prior; but since any prior on PX that is independent of the prior on ((qτj,Ij)j=1,…,s,σ,p) would disappear upon marginalization of the posterior of (PX,(qτj,Ij)j=1,…,s,σ,p) relatively to PX, we drop it in the sequel. Thus, it suffices to choose a prior distribution for qτj,Ij, for j=1,…,s, σ, and p.(1)
For the quantiles, we distinguish between the linear and nonlinear cases.
(a)
Linear case.qτj(X)=βτj(0)+XTβτj, j∈{1,…,s}, where XT denotes the transpose of X and βτj∈Rd.Due to (5), one should note that distinct quantiles differ only by the intercept βτj(0), so that all quantiles of interest have the same slope (βτ1=…=βτs=βp); then we consider the unknown vector parameter βˇ=(βτ1(0),…,βτs(0),βp) of dimension s+d with Gaussian prior,(9)βˇ~Ns+d0,Σˇ0,withΣˇ0 being a positive definite square matrix of dimension s+d.(b)
Nonlinear case. We define the function h on [0,1]×Xd, by h(τ,x)=qτ(x). We put a Gaussian process prior on h, h~GP(0,k), with a zero-mean function and a covariance function k: [0,1]2×X2d→R; following Sangnier et al. [26], we choose k to be decomposable and its form is given by(10)kτi,x,τj,x′=kxx,x′exp⁡-cτj-τk2,wherekx(x,x′)=exp⁡-bx-x′2, b and c are positive hyperparameters, and ∥·∥ is the Euclidean norm on Rd. Here, k((τi,x),(τj,x′)), i≠j, encodes the relation between the conditional quantiles qτi(x) and qτj(x′). As explained by Sangnier et al. [26], if c→0, the quantile curves are parallel so they do not cross. If c→+∞, the quantiles are learned independently and then they may cross. Thus, the choice of c is important to control the crossing.(2)
σ~IG(a0,b0) and a0>0,b0>0, where IG denotes the inverse Gamma distribution with positive hyperparameters a0 and b0.(3)
p~Beta(α0,λ0) and α0,λ0>0.


## 2.1. Model

We consider the quantile regression model(1)Yi=qpXi+ϵi,i=1,…,n,with (ϵi)i~iidALD(μ=0,σ,p), where μ, σ, and p, respectively, are the location, the scale, and the asymmetry parameters. This leads to assuming that the response variable Yi, given the covariable Xi, follows an ALD distribution:(2)Yi∣qpXi,σ,p~ALDqpXi,σ,p∀i=1,…,n,whose probability density function (p.d.f.) is given by(3)fyi;qpxi,σ,p=p1-pσexp⁡-ρpyi-qpxiσ,where ρp(u)=u(p-1u<0) is the check loss function. The explanatory variables Xi’s are supposed to be iid distributed according to an arbitrary continuous distribution, PX, whose support is Xd, X⊂R, d≥1. In what follows, the function qp(.) and the parameters σ, p, and PX are supposed to be unknown.In all subsequent sections, any quantity in bold represents a vector; capital letters denote random variables or vectors, whereas lowercase letters denote observed values of the corresponding random vectors.


## 2.2. Bayesian Procedure

### 2.2.1. Likelihood

To infer simultaneouslys distinct quantiles with a fully Bayesian approach, say qτ1,…,qτs with distinct orders τ1,…,τs, of the conditional distribution of Y∣X, one needs to characterize the likelihood through all these s unknown quantiles.As explained below, this is done first by partitioning the whole sample intos subsamples, which requires a sufficient number of observations, second by using the relation between any two quantiles of the ALD distribution, and third by rewriting the whole likelihood through s terms, where the j-th term only depends on qτj.(1)
Consider the following partition of the whole sample throughs subsamples:(4)XIj,YIj=Xi,Yii∈Ij,∀j∈1,…,s,whereIj=(j-1)r+1,…,jr with Ij⋂Ik=∅, ∀j≠k, and j,k∈{1,…,s}; suppose that r=n/s is integer, without loss of generality.(2)
To characterize the likelihood through all quantiles of interest, we use the relation that links anyτ-th quantile to the p-th quantile of the underlying ALD distribution (see Yu and Zhang [37] and Rodrigues and Fan [31]):(5)qτX=qpX+σgτ,p,whereg(τ,p)=1/1-plog⁡τ/p10<τ≤p-1/plog⁡(1-τ/1-p)1p<τ<1.(3)
From (5), we rewrite the model given by (1) as follows:(6)YIj=qτj,Ij-σgτj,p+ϵIj,∀j=1,…,s,where(7)YIj=Yijij∈Ij,qτj,Ij=qτjXijij∈Ij,gτj,p=gτj,p1ij∈Ijij∈Ij,ϵIj=ϵijij∈Ij.Then, the likelihood associated with the model given by (1) is rewritten as the product of s likelihoods; each one only depends on one subsample: (8)Lqτ1,I1,…,qτs,Is,σ,p;xI1,…,xIs,yI1,…,yIs=∏j=1s∏ij∈Ijp1-pσexp⁡-ρpyij-qτjxij-σgτj,pσ.


### 2.2.2. Priors

In this section, we specify the prior distributions for all unknown quantities. It is worth mentioning thatPX is unknown even it is not our primary interest, and therefore it should require a prior; but since any prior on PX that is independent of the prior on ((qτj,Ij)j=1,…,s,σ,p) would disappear upon marginalization of the posterior of (PX,(qτj,Ij)j=1,…,s,σ,p) relatively to PX, we drop it in the sequel. Thus, it suffices to choose a prior distribution for qτj,Ij, for j=1,…,s, σ, and p.(1)
For the quantiles, we distinguish between the linear and nonlinear cases.
(a)
Linear case.qτj(X)=βτj(0)+XTβτj, j∈{1,…,s}, where XT denotes the transpose of X and βτj∈Rd.Due to (5), one should note that distinct quantiles differ only by the intercept βτj(0), so that all quantiles of interest have the same slope (βτ1=…=βτs=βp); then we consider the unknown vector parameter βˇ=(βτ1(0),…,βτs(0),βp) of dimension s+d with Gaussian prior,(9)βˇ~Ns+d0,Σˇ0,withΣˇ0 being a positive definite square matrix of dimension s+d.(b)
Nonlinear case. We define the function h on [0,1]×Xd, by h(τ,x)=qτ(x). We put a Gaussian process prior on h, h~GP(0,k), with a zero-mean function and a covariance function k: [0,1]2×X2d→R; following Sangnier et al. [26], we choose k to be decomposable and its form is given by(10)kτi,x,τj,x′=kxx,x′exp⁡-cτj-τk2,wherekx(x,x′)=exp⁡-bx-x′2, b and c are positive hyperparameters, and ∥·∥ is the Euclidean norm on Rd. Here, k((τi,x),(τj,x′)), i≠j, encodes the relation between the conditional quantiles qτi(x) and qτj(x′). As explained by Sangnier et al. [26], if c→0, the quantile curves are parallel so they do not cross. If c→+∞, the quantiles are learned independently and then they may cross. Thus, the choice of c is important to control the crossing.(2)
σ~IG(a0,b0) and a0>0,b0>0, where IG denotes the inverse Gamma distribution with positive hyperparameters a0 and b0.(3)
p~Beta(α0,λ0) and α0,λ0>0.


## 2.2.1. Likelihood

To infer simultaneouslys distinct quantiles with a fully Bayesian approach, say qτ1,…,qτs with distinct orders τ1,…,τs, of the conditional distribution of Y∣X, one needs to characterize the likelihood through all these s unknown quantiles.As explained below, this is done first by partitioning the whole sample intos subsamples, which requires a sufficient number of observations, second by using the relation between any two quantiles of the ALD distribution, and third by rewriting the whole likelihood through s terms, where the j-th term only depends on qτj.(1)
Consider the following partition of the whole sample throughs subsamples:(4)XIj,YIj=Xi,Yii∈Ij,∀j∈1,…,s,whereIj=(j-1)r+1,…,jr with Ij⋂Ik=∅, ∀j≠k, and j,k∈{1,…,s}; suppose that r=n/s is integer, without loss of generality.(2)
To characterize the likelihood through all quantiles of interest, we use the relation that links anyτ-th quantile to the p-th quantile of the underlying ALD distribution (see Yu and Zhang [37] and Rodrigues and Fan [31]):(5)qτX=qpX+σgτ,p,whereg(τ,p)=1/1-plog⁡τ/p10<τ≤p-1/plog⁡(1-τ/1-p)1p<τ<1.(3)
From (5), we rewrite the model given by (1) as follows:(6)YIj=qτj,Ij-σgτj,p+ϵIj,∀j=1,…,s,where(7)YIj=Yijij∈Ij,qτj,Ij=qτjXijij∈Ij,gτj,p=gτj,p1ij∈Ijij∈Ij,ϵIj=ϵijij∈Ij.Then, the likelihood associated with the model given by (1) is rewritten as the product of s likelihoods; each one only depends on one subsample: (8)Lqτ1,I1,…,qτs,Is,σ,p;xI1,…,xIs,yI1,…,yIs=∏j=1s∏ij∈Ijp1-pσexp⁡-ρpyij-qτjxij-σgτj,pσ.


## 2.2.2. Priors

In this section, we specify the prior distributions for all unknown quantities. It is worth mentioning thatPX is unknown even it is not our primary interest, and therefore it should require a prior; but since any prior on PX that is independent of the prior on ((qτj,Ij)j=1,…,s,σ,p) would disappear upon marginalization of the posterior of (PX,(qτj,Ij)j=1,…,s,σ,p) relatively to PX, we drop it in the sequel. Thus, it suffices to choose a prior distribution for qτj,Ij, for j=1,…,s, σ, and p.(1)
For the quantiles, we distinguish between the linear and nonlinear cases.
(a)
Linear case.qτj(X)=βτj(0)+XTβτj, j∈{1,…,s}, where XT denotes the transpose of X and βτj∈Rd.Due to (5), one should note that distinct quantiles differ only by the intercept βτj(0), so that all quantiles of interest have the same slope (βτ1=…=βτs=βp); then we consider the unknown vector parameter βˇ=(βτ1(0),…,βτs(0),βp) of dimension s+d with Gaussian prior,(9)βˇ~Ns+d0,Σˇ0,withΣˇ0 being a positive definite square matrix of dimension s+d.(b)
Nonlinear case. We define the function h on [0,1]×Xd, by h(τ,x)=qτ(x). We put a Gaussian process prior on h, h~GP(0,k), with a zero-mean function and a covariance function k: [0,1]2×X2d→R; following Sangnier et al. [26], we choose k to be decomposable and its form is given by(10)kτi,x,τj,x′=kxx,x′exp⁡-cτj-τk2,wherekx(x,x′)=exp⁡-bx-x′2, b and c are positive hyperparameters, and ∥·∥ is the Euclidean norm on Rd. Here, k((τi,x),(τj,x′)), i≠j, encodes the relation between the conditional quantiles qτi(x) and qτj(x′). As explained by Sangnier et al. [26], if c→0, the quantile curves are parallel so they do not cross. If c→+∞, the quantiles are learned independently and then they may cross. Thus, the choice of c is important to control the crossing.(2)
σ~IG(a0,b0) and a0>0,b0>0, where IG denotes the inverse Gamma distribution with positive hyperparameters a0 and b0.(3)
p~Beta(α0,λ0) and α0,λ0>0.


## 3. Computations

To compute the full conditional distributions, we shall make use of the location-scale mixture of the ALD distribution (see Kozumi and Kobayashi [14]).Letω be an exponential latent variable with parameter 1/σ, denoted by E(1/σ), and let Z be a standard normal variable such that ω and Z are independent. If ϵ has an ALD distribution, ϵ~ALD(0,σ,p); then, it can be represented as a mixture of normal variable:(11)ϵ=γpω+δpσωz,where γp=1-2p/p(1-p) and δp2=2/p(1-p).Exploiting this augmented data structure, the model defined by the system of equations given by (6) admits, conditionally on ωIj=(ωi)i∈Ij~iidE(1/σ), the following Gaussian representation:(12)YIj=qτj,Ij-σgτj,p+γpωIj+δpσωIjzIj,j=1,…,s,where zIj=(zi)i∈Ij~iidN(0,1).


### 3.1. Full Conditional Distributions

From the equations given in (12) and the priors defined in Section 2.2.2, we follow Bernardi et al. [15] to derive the full conditional distributions of ωn=(ω1,…,ωn) and (qτj,Ij)j=1,…,s. Yet, the full conditionals of σ and p are not tractable and this is carried out by including a Metropolis-Hasting step to the Gibbs sampler algorithm.(1)
For allj=1,…,s, as well as i∈Ij, set νi=ωi-1, and considering the distribution ωi~E(1/σ) as a prior on ωi,(13)νi∣yi,qτjxi,σ,p~IGaussΨi,ϕpwith(14)Ψi=γp2+2δp2yi-qτjxi+σgτj,p2,ϕp=γp2+2δp2δp2σ,whereIGauss stands for the Inverse Gaussian distribution with Ψi>0 and ϕp>0 as location and shape parameters.(2)
The conjugate Gaussian prior on quantiles provides Gaussian full conditional distributions in both linear and nonlinear cases.
(a)
Linear case. Let us introduce some notations:
(i)
Yn=(Yi)1≤i≤n,(ii)
Xˇ=(Xˇi,l)1≤i≤n,1≤l≤s+d is the n×(s+d) design matrix defined by(15)Xˇi,l=1i∈Ilifl∈1,…s,Xi,l-sifl∈s+1,…,s+d,(iii)
gn=(g(τj,p)1i∈Ij)1≤j≤s.This allows to rewrite the system of equations (12) into the following format:(16)Yn=Xˇβˇ-σgn+γpωn+ϵn∗,withϵn∗~Nn(0,Σ∗) and where Σ∗=diag(δp2σωn). With the zero-mean Gaussian prior distribution, the full conditional distribution on βˇ is then Gaussian:(17)βˇ∣xˇ,yn,ωn,σ,p~Ns+dμ^βˇ,Σ^βˇwith(18)μ^βˇ=xˇTΣ∗xˇ+Σˇ0-1-1xˇTΣ∗yn+σgn-γpωn,Σ^βˇ=xˇTΣ∗xˇ+Σˇ0-1-1.(b)
Nonlinear case. We shall use other extra notations: (19)hn=hτj,Ijj=1,…,s=hτj,xii∈Ijj=1,…,s,knτ,x=kτj,xijij∈Ijj=1,…,s.This allows to rewrite the system of equations (12) into the following vector format:(20)Yn=hn-σgn+γpωn+ϵn∗.Combining (20) with the Gaussian process prior on h leads to the following joint distribution of (Yn,hn) conditional on (x,ωn,σ,p):(21)Ynhn~N2n-σgn+γpωn0n,knτ,x+Σ∗knτ,xknτ,xknτ,x.Finally, classical calculations lead to the desired full conditional distribution:(22)hn∣x,yn,ωn,σ,p~Nμ^h,Σ^h,where(23)μ^h=knτ,xknτ,x+Σ∗-1yn-γpωn+σgn,Σ^h=knτ,x-knτ,xknτ,x+Σ∗-1knτ,x.(3)
The full conditional distribution,πσ of σ∣x,ωn,p, is proportional to(24)1σn+a0+1exp⁡-12Yn-qτ,n+σgn-γpωnTΣ∗-1Yn-qτ,n+σgn-γpωn,whereqτ,n=qτj,Ijj=1,…,s.(4)
The full conditional,πp of p∣x,ωn,σ, is proportional to(25)pα01-pλ0exp⁡-12Yn-qτ,n+σgn-γpωnTΣ∗-1Yn-qτ,n+σgn-γpωn.


### 3.2. Algorithm

Due to (24) and (25), it is not possible to generate σ and p directly from their full conditional distribution. Therefore, they are simulated by incorporating a random walk Metropolis-Hastings step within Gibbs, as described in Algorithm 1.Algorithm 1:Metropolis-Hastings step.
1:Initializationt=0: σ(0)=σ0, p(0)=p0. 2:Loopt=t+1: 1.    (i)σprop|σt~N]0,∞]σt,sdσ, (ii)ασ=min⁡1,PσQσ (iii)u1~U[0,1]⇒σ(t+1))=σpropifα>u1σ(t)otherwise 2.   (i)pprop|p(t)~N]0,1[(p(t),sdp), (ii)αp=min⁡1,PpQp (iii)u2~U[0,1]⇒p(t+1)=ppropifα>u2p(t)otherwiseFor a subsetE⊂R, denote by NE(·,·) the truncated version on E of the corresponding Gaussian distribution. Note that the choice of the truncated Gaussian distribution is classical.It is worth to mention that the values of the scale parameterssdσ and sdp are calibrated to achieve the equilibrium of the random walk Metropolis-Hastings step quickly; in fact, they are chosen neither too small nor too large so that the acceptance rate becomes practically stable.The conditional posterior ratios ofσ and p are given by(26)Pσ=πσσprop∣yn,xn,qτj,Ijj=1,…,st+1,ωnt+1,ptπσσt∣yn,xn,qτj,Ijj=1,…,st+1,ωnt+1,ptPp=πppprop∣yn,xn,qτj,Ijj=1,…,st+1,ωnt+1,ptπppt∣yn,xn,qτj,Ijj=1,…,st+1,ωnt+1,pt,where πσ(·∣·) and πp(·∣·), are respectively, given up to a constant in (24) and (25), and the transition probabilities are given by(27)Qσ=fN0,∞σt,sdσfN0,∞σprop,sdσ=1-ϕ-σprop/sdσ1-ϕ-σt/sdσ,Qp=fN0,1pt,sdpfN0,1pprop,sdp=ϕ1-pprop/sdp-ϕ-pprop/sdpϕ1-pt/sdp-ϕ-pt/sdp,where fNE(·,·) denotes the pdf of NE(·,·) and ϕ(.) denotes the cumulative distribution function of the standard normal distribution.


## 3.1. Full Conditional Distributions

From the equations given in (12) and the priors defined in Section 2.2.2, we follow Bernardi et al. [15] to derive the full conditional distributions of ωn=(ω1,…,ωn) and (qτj,Ij)j=1,…,s. Yet, the full conditionals of σ and p are not tractable and this is carried out by including a Metropolis-Hasting step to the Gibbs sampler algorithm.(1)
For allj=1,…,s, as well as i∈Ij, set νi=ωi-1, and considering the distribution ωi~E(1/σ) as a prior on ωi,(13)νi∣yi,qτjxi,σ,p~IGaussΨi,ϕpwith(14)Ψi=γp2+2δp2yi-qτjxi+σgτj,p2,ϕp=γp2+2δp2δp2σ,whereIGauss stands for the Inverse Gaussian distribution with Ψi>0 and ϕp>0 as location and shape parameters.(2)
The conjugate Gaussian prior on quantiles provides Gaussian full conditional distributions in both linear and nonlinear cases.
(a)
Linear case. Let us introduce some notations:
(i)
Yn=(Yi)1≤i≤n,(ii)
Xˇ=(Xˇi,l)1≤i≤n,1≤l≤s+d is the n×(s+d) design matrix defined by(15)Xˇi,l=1i∈Ilifl∈1,…s,Xi,l-sifl∈s+1,…,s+d,(iii)
gn=(g(τj,p)1i∈Ij)1≤j≤s.This allows to rewrite the system of equations (12) into the following format:(16)Yn=Xˇβˇ-σgn+γpωn+ϵn∗,withϵn∗~Nn(0,Σ∗) and where Σ∗=diag(δp2σωn). With the zero-mean Gaussian prior distribution, the full conditional distribution on βˇ is then Gaussian:(17)βˇ∣xˇ,yn,ωn,σ,p~Ns+dμ^βˇ,Σ^βˇwith(18)μ^βˇ=xˇTΣ∗xˇ+Σˇ0-1-1xˇTΣ∗yn+σgn-γpωn,Σ^βˇ=xˇTΣ∗xˇ+Σˇ0-1-1.(b)
Nonlinear case. We shall use other extra notations: (19)hn=hτj,Ijj=1,…,s=hτj,xii∈Ijj=1,…,s,knτ,x=kτj,xijij∈Ijj=1,…,s.This allows to rewrite the system of equations (12) into the following vector format:(20)Yn=hn-σgn+γpωn+ϵn∗.Combining (20) with the Gaussian process prior on h leads to the following joint distribution of (Yn,hn) conditional on (x,ωn,σ,p):(21)Ynhn~N2n-σgn+γpωn0n,knτ,x+Σ∗knτ,xknτ,xknτ,x.Finally, classical calculations lead to the desired full conditional distribution:(22)hn∣x,yn,ωn,σ,p~Nμ^h,Σ^h,where(23)μ^h=knτ,xknτ,x+Σ∗-1yn-γpωn+σgn,Σ^h=knτ,x-knτ,xknτ,x+Σ∗-1knτ,x.(3)
The full conditional distribution,πσ of σ∣x,ωn,p, is proportional to(24)1σn+a0+1exp⁡-12Yn-qτ,n+σgn-γpωnTΣ∗-1Yn-qτ,n+σgn-γpωn,whereqτ,n=qτj,Ijj=1,…,s.(4)
The full conditional,πp of p∣x,ωn,σ, is proportional to(25)pα01-pλ0exp⁡-12Yn-qτ,n+σgn-γpωnTΣ∗-1Yn-qτ,n+σgn-γpωn.


## 3.2. Algorithm

Due to (24) and (25), it is not possible to generate σ and p directly from their full conditional distribution. Therefore, they are simulated by incorporating a random walk Metropolis-Hastings step within Gibbs, as described in Algorithm 1.Algorithm 1:Metropolis-Hastings step.
1:Initializationt=0: σ(0)=σ0, p(0)=p0. 2:Loopt=t+1: 1.    (i)σprop|σt~N]0,∞]σt,sdσ, (ii)ασ=min⁡1,PσQσ (iii)u1~U[0,1]⇒σ(t+1))=σpropifα>u1σ(t)otherwise 2.   (i)pprop|p(t)~N]0,1[(p(t),sdp), (ii)αp=min⁡1,PpQp (iii)u2~U[0,1]⇒p(t+1)=ppropifα>u2p(t)otherwiseFor a subsetE⊂R, denote by NE(·,·) the truncated version on E of the corresponding Gaussian distribution. Note that the choice of the truncated Gaussian distribution is classical.It is worth to mention that the values of the scale parameterssdσ and sdp are calibrated to achieve the equilibrium of the random walk Metropolis-Hastings step quickly; in fact, they are chosen neither too small nor too large so that the acceptance rate becomes practically stable.The conditional posterior ratios ofσ and p are given by(26)Pσ=πσσprop∣yn,xn,qτj,Ijj=1,…,st+1,ωnt+1,ptπσσt∣yn,xn,qτj,Ijj=1,…,st+1,ωnt+1,ptPp=πppprop∣yn,xn,qτj,Ijj=1,…,st+1,ωnt+1,ptπppt∣yn,xn,qτj,Ijj=1,…,st+1,ωnt+1,pt,where πσ(·∣·) and πp(·∣·), are respectively, given up to a constant in (24) and (25), and the transition probabilities are given by(27)Qσ=fN0,∞σt,sdσfN0,∞σprop,sdσ=1-ϕ-σprop/sdσ1-ϕ-σt/sdσ,Qp=fN0,1pt,sdpfN0,1pprop,sdp=ϕ1-pprop/sdp-ϕ-pprop/sdpϕ1-pt/sdp-ϕ-pt/sdp,where fNE(·,·) denotes the pdf of NE(·,·) and ϕ(.) denotes the cumulative distribution function of the standard normal distribution.


## 4. Simulation Study

In this section, we study the performance of our method in both linear and nonlinear quantile regressions. For the model given by (1), we shall consider three different designs for the p-th quantile:(1)
Univariate linear quantile:qp(X)=1+2X, with X~U[-1,1].(2)
Multivariate linear quantile:qp(X)=1+XTβp, with βp∈Rd, d=10, and either Xl~U[-1,1] or Xl~N(0,1), ∀l∈{1,…,d}.(3)
Nonlinear quantile: forX~U[0,1], (28)qpX=cos⁡52πXexp⁡-32X+14exp⁡2X-0.5-exp⁡-114+X21X<0.5+14exp⁡-2X-0.5-14+X2exp⁡-1-12-X1X≥0.5.For all designs, we generate independently 300 observations issued from the model defined in (1). All runs of Metropolis-Hastings within Gibbs algorithm consist in 20000 iterations, one-third of which are burn-in. The prior hyperparameters are chosen as follows: for the inverse Gamma of σ, a0=1 and b0=0.01; for the Beta prior of p, α0=2 and λ0=2 and Σˇ0 is set to be the identity matrix for linear quantile and c=0.1 and b=5 for the nonlinear case. The choice of c=0.1 and b=5 is typical; indeed, whatever the value of τ is, c=0.1 minimizes the empirical root mean integrated square error, RMISE; that is, c=argmin1/n∑i=1n(qτ(Xi)-q^τ(Xi))2, where q^τ stands for the posterior mean quantile regression. In order to test the robustness of our procedure with respect to the model parameters, different values of σ and p are considered for the three designs.The first design is a straightforward example and is carried out just to check the convergence of the algorithm from different tools:R^ of Gelman and Rubin diagnostic, the autocorrelation analysis, and the posterior plots of the various parameters.Through the second design, we use the crossing loss criterion (see Sangnier et al. [26]) to compare the performance, against crossing, of our method, denoted by “SBQR" with other approaches: the frequentist single quantile method of Koenker and Bassett [1], denoted by “K&B”, the Bayesian single quantile regression method of Yu and Moyeed [13], denoted by “Y&M", and the simultaneous Bayesian method of Reich and Smith [38], denoted by “R&S”; in addition, we use the RMISE to evaluate the estimation performance. These other methods are performed using available codes in R Core Team (2017):rq function available inquantreg package Koenker [39] for “K&B”,bayesQR function in bayesQR package Benoit et al. [40] for “Y&M", andqreg function inBSquare package Smith et al. [41] for “R&S".For design 3, we compare our “SBQR" method with both the nonparametric quantile regression method of Muggeo et al. [32], denoted by “M&ST", and the simultaneous noncrossing method of Rodrigues and Fan [31], denoted by “F&R”. We have implemented “M&ST" with quantregGrowth R package (see Muggeo [42]) and “F&R” with an own made code in R. We use features of the RMISE criterion to show how well “SBQR” performs among the other considered methods.We note that the number of observations,n=300, is chosen carefully to illustrate the intended objectives in all designs.


### 4.1. Design 1: Univariate Linear Quantile Regression

The iid sampleϵ1,…,ϵn is generated according to ALD(0,σ=0.1,p=0.25). We propose to infer three quantiles that are close, namely, the quantiles of order τ=0.2,0.3 and 0.4, respectively.We fixsdσ=sdp=0.1. To evaluate the convergence of our algorithm, we use three different seeds and parameter starting values to run three different chains and calculate R^ of Gelman convergence diagnostic. Besides, we use other convergence diagnostics such as the autocorrelation analysis and the posterior plots.As shown in the top panel of Figures1, 2, and 3, all posterior distributions shrink at the true parameters value. Furthermore, in the middle panel of Figures 1, 2, and 3, the decrease of the empirical autocorrelation of posterior samples proves that the underlying chains are stationary. The bottom panels of Figures 1, 2, and 3 show that R^ goes to 1 through the iterations, which confirms the convergence of the algorithm.Figure 1
Trace and density plot (top), autocorrelation plot (middle), andR^ evolution through iterations (bottom) of σ^ (left panel) and p^ (right panel).Figure 2
Trace and density plot (top), autocorrelation plot (middle), andR^ evolution (bottom) of the intercept β^0.2(0) (left panel) and the slope β^0.3(0) of the 0.2-th and 0.3-th quantile, respectively.Figure 3
Trace and density plot (top), autocorrelation plot (middle), andR^ evolution (bottom) of the intercept β^0.4(0) (left panel) of the 0.4-th quantile and the slope β^p.


### 4.2. Design 2: Multivariate Linear Quantile Regression

The second design is dedicated to the multivariate linear case; hence, we consider the model given by (1) with ϵ1,…,ϵn~iidALD(μ=0,σ=0.5,p=0.25), X∈Rd, d=10, Xk~U[-1,1],k=1,…,d, and βp=(1.6, 2.2, 2.8, 3.4, 4, 4.5, 5.1, 5.7, 6.3, 6.9).As commonly known, crossing quantiles is a practical problem that often occurs when there is a large number of covariates. We propose to infer two close quantiles of orderτ1=0.1 and τ2=0.12 and to study the crossing throughout the following four methods: “K&B", “Y&M", “R&S", and “SBQR".To achieve the desired posterior distribution through MCMC methods, we perform Y&M and R&S with different number of iterations: 1000 for Y&M and 10000 for R&S. For “R&S”, we use the logistic base distribution with 4 basis functions. For “SBQR", we fixsdσ=sdp=0.01 and Σˇ0 to be the identity matrix.To compare the methods, we make use of the crossing loss criterion (see Sangnier et al. [26]) that measures how far q^0.12(X) goes below q^0.10(X):(29)crossingloss=1n∑i=1nmax⁡0,q^0.10Xi-q^0.12Xi.For a given approach, the less the crossing loss is, the better is the method. As shown in Table1, the crossing loss, by “K&B", is significantly of 0.43%, which corresponds to 34 data points of q^0.1(X) which are above q^0.12(X). This percentage is considerably weakened when applying the separate Bayesian method “Y&M" (0.38%) but still has crossing quantiles (26 data points of q^0.1(X) are above q^0.12(X)). However, for simultaneous estimation methods, as our proposed “SBQR" or “R&S” methods, the crossing loss becomes zero; this means that the simultaneous quantile estimation has the potential to make quantile crossing vanish. Moreover, the estimated quantiles are in the right order. While simultaneous approaches control the monotonicity property of quantiles in a certain sense, separate approaches do not, and they provide more than 8% of violation according to results previously discussed (11.33% of violation by “K&B" and 8.66% of violation by “Y&M").Table 1
Table of criteria: crossing loss and RMISE.
Method nb. crossing crossing loss RMISE (0.10) RMISE (0.12) “K&B" 34 0.0043 0.3565 0.3218 “Y&M" 26 0.0038 0.3344 0.3362 “R&S" 0 0 0.2449 0.2849 “SBQR" 0 0 0.1989 0.2383To go further, we compute the RMISE for each quantile order and for all methods. The “SBQR" has roughly the smallest RMISE in both quantile levelsτ=0.1 and τ=0.12. The fact is that “R&S” may oversmooth when estimating the quantiles simultaneously and here, for this linear case, the oversmoothness impacts the flexibility of the method (high RMISE).To see if the covariable support has an impact on results, we consider another simulation set in whichX~Nd(0d,Idd), that is, support(X)=Rd. Table 2 shows that “SBQR" behaves like in the previous example: a zero crossing loss and the smallest RMISE among all the methods.Table 2
Table of criteria: crossing loss and RMISE for normal covariate case.
Method nb. crossing crossing loss RMISE (0.10) RMISE (0.12) “K&B" 38 0.0037 0.2356 0.2151 “Y&M" 5 0.0005 0.4569 0.4391 “R&S" 0 0 0.2865 0.3189 “SBQR" 0 0 0.2218 0.2102We also consider another case with a different pair of quantile orders:τ1=0.7 and τ2=0.8. We turn back to the support [-1,1]d for X. Table 3 shows similar results as the ones obtained for τ1=0.10 and τ2=0.12. Thus, “SBQR” still has the best behavior among the other methods in terms of crossing loss and RMISE.Table 3
Table of criteria: crossing loss and RMISE for uniform covariate case.
Method nb. crossing crossing loss RMISE (0.7) RMISE (0.8) “K&B" 6 0.0041 0.4255 0.7177 “Y&M" 1 0.0009 0.403 0.7013 “R&S" 0 0 0.224 0.3791 “SBQR" 0 0 0.2065 0.2708It should be noted that the estimation ofσ and p by “SBQR" is quite good, since their estimated values are near the true ones in the different treated cases.


### 4.3. Design 3: Nonparametric Quantile Regression

Considering the third design withX~U[0,1] and ϵ~ALD(μ=0,σ=0.05,p=0.75), we are interested in estimating quantile functions for orders τ=0.10, 0.12, 0.15, and 0.20. We fix sdp=0.05 and sdσ=0.005 and we compare our “SBQR" method with two others: the “M&ST” method with three-order cubic B-splines and the “F&R" approach with smoothness parameter value equal to 0.1.For each value ofτ, we evaluate the performance of these methods through the RMISE criterion. Table 4 shows that the RMISE values are significantly smaller for “SBQR" than the ones for “M&ST” and “F&R”. It is worth noting, therefore, that our method is significantly better than the two others in quantiles estimation; besides, it provides a reasonable estimation of σ (σ^=0.0535) and p (p^=0.7744).Table 4
RMISE at different quantile levels computed for M&ST, F&R, and SBQR methods.
Method RMISE (0.10) RMISE (0.12) RMISE (0.15) RMISE (0.20) “M&ST" 0.1416 0.1397 0.1072 0.0774 “F&R" 0.1520 0.1033 0.0821 0.0664 “SBQR" 0.0692 0.0768 0.0777 0.0660Figure4 gives the quantile curves estimators. While our “SBQR" method (left panel) provides quantile curves estimates that are close to the true ones (red-dashed lines), the 0.1-th quantile curve estimate given by “F&R" method (middle panel) is fairly distant from the true curve, especially when x∈[0.2,0.6]. The same happens for “M&ST" method when x∈[0.7,1]. However, there is no scarred crossing by any of these three methods, since they are tackling simultaneous quantile estimation techniques.Figure 4
Estimated quantile curves (black solid lines) against the true ones (red dashed lines) for SBQR (left panel), F&R (middle panel), and M&ST (right panel) methods.The issue that makes “F&R" method less flexible in simultaneous quantiles fitting is that, in the second stage, the final estimators are quite affected by the first stage output that may be badly estimated. For “M&ST", the estimators are constructed iteratively, when solving the minimization problem, by adding constraints so that each subsequent quantile function does not cross with the previous one; this may cause overestimation of quantile curves. Our “SBQR" method does not suffer from these mentioned disadvantages, which explains its significant flexibility.


## 4.1. Design 1: Univariate Linear Quantile Regression

The iid sampleϵ1,…,ϵn is generated according to ALD(0,σ=0.1,p=0.25). We propose to infer three quantiles that are close, namely, the quantiles of order τ=0.2,0.3 and 0.4, respectively.We fixsdσ=sdp=0.1. To evaluate the convergence of our algorithm, we use three different seeds and parameter starting values to run three different chains and calculate R^ of Gelman convergence diagnostic. Besides, we use other convergence diagnostics such as the autocorrelation analysis and the posterior plots.As shown in the top panel of Figures1, 2, and 3, all posterior distributions shrink at the true parameters value. Furthermore, in the middle panel of Figures 1, 2, and 3, the decrease of the empirical autocorrelation of posterior samples proves that the underlying chains are stationary. The bottom panels of Figures 1, 2, and 3 show that R^ goes to 1 through the iterations, which confirms the convergence of the algorithm.Figure 1
Trace and density plot (top), autocorrelation plot (middle), andR^ evolution through iterations (bottom) of σ^ (left panel) and p^ (right panel).Figure 2
Trace and density plot (top), autocorrelation plot (middle), andR^ evolution (bottom) of the intercept β^0.2(0) (left panel) and the slope β^0.3(0) of the 0.2-th and 0.3-th quantile, respectively.Figure 3
Trace and density plot (top), autocorrelation plot (middle), andR^ evolution (bottom) of the intercept β^0.4(0) (left panel) of the 0.4-th quantile and the slope β^p.


## 4.2. Design 2: Multivariate Linear Quantile Regression

The second design is dedicated to the multivariate linear case; hence, we consider the model given by (1) with ϵ1,…,ϵn~iidALD(μ=0,σ=0.5,p=0.25), X∈Rd, d=10, Xk~U[-1,1],k=1,…,d, and βp=(1.6, 2.2, 2.8, 3.4, 4, 4.5, 5.1, 5.7, 6.3, 6.9).As commonly known, crossing quantiles is a practical problem that often occurs when there is a large number of covariates. We propose to infer two close quantiles of orderτ1=0.1 and τ2=0.12 and to study the crossing throughout the following four methods: “K&B", “Y&M", “R&S", and “SBQR".To achieve the desired posterior distribution through MCMC methods, we perform Y&M and R&S with different number of iterations: 1000 for Y&M and 10000 for R&S. For “R&S”, we use the logistic base distribution with 4 basis functions. For “SBQR", we fixsdσ=sdp=0.01 and Σˇ0 to be the identity matrix.To compare the methods, we make use of the crossing loss criterion (see Sangnier et al. [26]) that measures how far q^0.12(X) goes below q^0.10(X):(29)crossingloss=1n∑i=1nmax⁡0,q^0.10Xi-q^0.12Xi.For a given approach, the less the crossing loss is, the better is the method. As shown in Table1, the crossing loss, by “K&B", is significantly of 0.43%, which corresponds to 34 data points of q^0.1(X) which are above q^0.12(X). This percentage is considerably weakened when applying the separate Bayesian method “Y&M" (0.38%) but still has crossing quantiles (26 data points of q^0.1(X) are above q^0.12(X)). However, for simultaneous estimation methods, as our proposed “SBQR" or “R&S” methods, the crossing loss becomes zero; this means that the simultaneous quantile estimation has the potential to make quantile crossing vanish. Moreover, the estimated quantiles are in the right order. While simultaneous approaches control the monotonicity property of quantiles in a certain sense, separate approaches do not, and they provide more than 8% of violation according to results previously discussed (11.33% of violation by “K&B" and 8.66% of violation by “Y&M").Table 1
Table of criteria: crossing loss and RMISE.
Method nb. crossing crossing loss RMISE (0.10) RMISE (0.12) “K&B" 34 0.0043 0.3565 0.3218 “Y&M" 26 0.0038 0.3344 0.3362 “R&S" 0 0 0.2449 0.2849 “SBQR" 0 0 0.1989 0.2383To go further, we compute the RMISE for each quantile order and for all methods. The “SBQR" has roughly the smallest RMISE in both quantile levelsτ=0.1 and τ=0.12. The fact is that “R&S” may oversmooth when estimating the quantiles simultaneously and here, for this linear case, the oversmoothness impacts the flexibility of the method (high RMISE).To see if the covariable support has an impact on results, we consider another simulation set in whichX~Nd(0d,Idd), that is, support(X)=Rd. Table 2 shows that “SBQR" behaves like in the previous example: a zero crossing loss and the smallest RMISE among all the methods.Table 2
Table of criteria: crossing loss and RMISE for normal covariate case.
Method nb. crossing crossing loss RMISE (0.10) RMISE (0.12) “K&B" 38 0.0037 0.2356 0.2151 “Y&M" 5 0.0005 0.4569 0.4391 “R&S" 0 0 0.2865 0.3189 “SBQR" 0 0 0.2218 0.2102We also consider another case with a different pair of quantile orders:τ1=0.7 and τ2=0.8. We turn back to the support [-1,1]d for X. Table 3 shows similar results as the ones obtained for τ1=0.10 and τ2=0.12. Thus, “SBQR” still has the best behavior among the other methods in terms of crossing loss and RMISE.Table 3
Table of criteria: crossing loss and RMISE for uniform covariate case.
Method nb. crossing crossing loss RMISE (0.7) RMISE (0.8) “K&B" 6 0.0041 0.4255 0.7177 “Y&M" 1 0.0009 0.403 0.7013 “R&S" 0 0 0.224 0.3791 “SBQR" 0 0 0.2065 0.2708It should be noted that the estimation ofσ and p by “SBQR" is quite good, since their estimated values are near the true ones in the different treated cases.


## 4.3. Design 3: Nonparametric Quantile Regression

Considering the third design withX~U[0,1] and ϵ~ALD(μ=0,σ=0.05,p=0.75), we are interested in estimating quantile functions for orders τ=0.10, 0.12, 0.15, and 0.20. We fix sdp=0.05 and sdσ=0.005 and we compare our “SBQR" method with two others: the “M&ST” method with three-order cubic B-splines and the “F&R" approach with smoothness parameter value equal to 0.1.For each value ofτ, we evaluate the performance of these methods through the RMISE criterion. Table 4 shows that the RMISE values are significantly smaller for “SBQR" than the ones for “M&ST” and “F&R”. It is worth noting, therefore, that our method is significantly better than the two others in quantiles estimation; besides, it provides a reasonable estimation of σ (σ^=0.0535) and p (p^=0.7744).Table 4
RMISE at different quantile levels computed for M&ST, F&R, and SBQR methods.
Method RMISE (0.10) RMISE (0.12) RMISE (0.15) RMISE (0.20) “M&ST" 0.1416 0.1397 0.1072 0.0774 “F&R" 0.1520 0.1033 0.0821 0.0664 “SBQR" 0.0692 0.0768 0.0777 0.0660Figure4 gives the quantile curves estimators. While our “SBQR" method (left panel) provides quantile curves estimates that are close to the true ones (red-dashed lines), the 0.1-th quantile curve estimate given by “F&R" method (middle panel) is fairly distant from the true curve, especially when x∈[0.2,0.6]. The same happens for “M&ST" method when x∈[0.7,1]. However, there is no scarred crossing by any of these three methods, since they are tackling simultaneous quantile estimation techniques.Figure 4
Estimated quantile curves (black solid lines) against the true ones (red dashed lines) for SBQR (left panel), F&R (middle panel), and M&ST (right panel) methods.The issue that makes “F&R" method less flexible in simultaneous quantiles fitting is that, in the second stage, the final estimators are quite affected by the first stage output that may be badly estimated. For “M&ST", the estimators are constructed iteratively, when solving the minimization problem, by adding constraints so that each subsequent quantile function does not cross with the previous one; this may cause overestimation of quantile curves. Our “SBQR" method does not suffer from these mentioned disadvantages, which explains its significant flexibility.


## 5. Conclusion

Our proposed estimation procedure, “SBQR", for simultaneous Bayesian quantile regression guarantees the fundamental property of noncrossing. Assuming that the ALD is the underlying data distribution, this method enables characterizing the likelihood function by all quantiles of interest using the relation between two distinct quantiles. Using the location-scale normal mixture representation of the ALD distribution, we develop a Metropolis-Hastings within Gibbs algorithm to implement our method. Our simulation studies show good results that reflect the good performance of the method and the convergence of the algorithm. Against the crossing problem of estimated quantiles, our method has good performance compared with single quantile estimation methods like Koenker and Bassett [1] and Yu and Moyeed [13] methods. From the RMISE point of view, “SBQR" is very competitive in both linear and nonlinear quantile regression cases.As future perspectives, we envisage extending this work to deal with any conditional distribution ofY∣X.



---

*Source: 8610723-2019-06-02.xml*