A New Family of Error Distributions for Bayesian Quantile Regression 7 1 YifeiYanandAthanasiosKottas∗ 0 2 b e F 9 Abstract ] E We propose a new family of error distributions for model-based quantile regression, which is M constructed through a structured mixture of normal distributions. The construction enables fix- . t a ing specific percentiles of the distribution while, at the same time, allowing for varying mode, t s skewness and tail behavior. It thus overcomes the severe limitation of the asymmetric Laplace [ 2 distribution–themostcommonlyusederrormodelforparametricquantileregression–forwhich v 6 theskewnessoftheerrordensityisfullyspecifiedwhenaparticularpercentileisfixed. Wede- 6 6 velopaBayesianformulationfortheproposedquantileregressionmodel,includingconditional 5 lassoregularizedquantileregressionbasedonahierarchicalLaplacepriorfortheregressionco- 0 . 1 efficients,andaTobitquantileregressionmodel. PosteriorinferenceisimplementedviaMarkov 0 7 ChainMonteCarlomethods. TheflexibilityofthenewmodelrelativetotheasymmetricLaplace 1 : distributionisstudiedthroughrelevantmodelproperties,andthroughasimulationexperimentto v i X comparethetwoerrordistributionsinregularizedquantileregression. Moreover,modelperfor- r manceinlinearquantileregression,regularizedquantileregression,andTobitquantileregression a isillustratedwithdataexamplesthathavebeenpreviouslyconsideredintheliterature. Keywords. AsymmetricLaplacedistribution;MarkovchainMonteCarlo;Regularizedquantile regression;Skewnormaldistribution;Tobitquantileregression. ∗Y.Yan(yifeiyan@soe.ucsc.edu)isPh.D.student,andA.Kottas(thanos@soe.ucsc.edu)isProfessorofStatistics, DepartmentofAppliedMathematicsandStatistics,UniversityofCalifornia,SantaCruz,CA,95064,USA. 1 1 Introduction Quantileregressionoffersapracticallyimportantalternativetotraditionalmeanregression,andforms an area with a rapidly increasing literature. Parametric quantile regression models are almost exclu- sivelybuiltfromtheasymmetricLaplace(AL)distributionthedensityofwhichis (cid:26) (cid:27) p(1−p) 1 fAL(y | µ,σ) = exp − ρ (y−µ) , y ∈ R (1) p σ σ p where ρ (u) = u[p − I(u < 0)], with I(·) denoting the indicator function. Here, σ > 0 is a p scale parameter, p ∈ (0,1), and µ ∈ R corresponds to the pth percentile, (cid:82)µ fAL(y | µ,σ)dy = −∞ p p. Hence, a model for pth quantile regression can be developed by expressing µ as a function of available covariates x, for instance, µ = xTβ yields a linear quantile regression structure. Note that maximizing the likelihood with respect to β under an AL response distribution corresponds to minimizing for β the check loss function, (cid:80)n ρ (y − xTβ), used for classical semiparametric i=1 p i i estimationinlinearquantileregression(Koenker,2005). The AL distribution is receiving increasing attention in the Bayesian literature, originating from workoninferenceforlinearquantileregression(YuandMoyeed,2001;Tsionas,2003). Particularly relevanttotheBayesianframeworkarethedifferentmixturerepresentationsofthedistribution(Kotz et al., 2001), which have been exploited to construct posterior simulation algorithms (Kozumi and Kobayashi, 2011), as well as to explore different modeling scenarios; see, for instance, Lum and Gelfand(2012)andWaldmannetal.(2013). However,theALdistributionhassubstantiallimitationsasanerrormodelforquantileregression. Most striking is that the skewness of the error density is fully determined when a specific percentile ischosen, thatis, whenpisfixed. Inparticular, theerrordensityissymmetricinthecaseofmedian regression,sinceforp = 0.5,theALreducestotheLaplacedistribution. Moreover,themodeofthe errordistributionisatzero,foranyp,whichresultsinrigiderrordensitytailsforextremepercentiles. The literature includes Bayesian nonparametric models for the error distribution in the special case of median regression (Walker and Mallick, 1999; Kottas and Gelfand, 2001; Hanson and John- son, 2002) and in general quantile regression (Kottas and Krnjajic´, 2009; Reich et al., 2010). The 2 Bayes nonparametrics literature has also explored inference methods for simultaneous quantile regression (Taddy and Kottas, 2010; Tokdar and Kadane, 2012; Reich and Smith, 2013). However, work on parametric alternatives to AL quantile regression errors is limited, and the existing models do not overcome all the limitations discussed above. For instance, although the class of skew distributions studied in Wichitaksorn et al. (2014) includes the AL as a special case, it shares the same restriction with the AL as a quantile regression error model in that it has a single parameter that controls both skewness and percentiles. Zhu and Zinde-Walsh (2009) and Zhu and Galbraith (2011) explored the family of asymmetric exponential power distributions, which does not include the AL distribution. For a fixed probability p, the density function has four free parameters and allows for different decay rates in the left and the right tails. However, similar to the AL, the mode of the distributionisfixedatthequantileµbyconstruction. Moreflexibleparametricquantileregressionerrormodelsarearguablyusefulbothtoexpandthe inferential scope of the asymmetric Laplace in the standard quantile regression setting, as well as to provide building blocks for model development under more complex data structures. The limited scopeofresultsinthisdirectionmaybeattributedtothechallengeofdefiningsufficientlyflexibledis- tributionsthatareparameterizedbypercentilesand,atthesametime,allowforpracticablemodeling andinferencemethods. Seekingtofillthisgap,weproposeanewfamilyofdistributionsthatisparameterizedintermsof percentiles,andovercomestherestrictiveaspectsoftheALdistribution. Thedistributionisdeveloped constructively through an extension of an AL mixture representation. In particular, we introduce a shape parameter to obtain a distribution that has more flexible skewness and tail behaviour than the ALdistribution,whileretainingitasaspecialcaseofthenewmodel. Thelatterenablesconnections with the check loss function which are useful in studying the utility of the new model in the context of regularized quantile regression. Owing to its hierarchical mixture representation, the proposed distribution preserves the important feature of ready to implement posterior inference for Bayesian quantileregression. In Section 2, we develop the new distribution and discuss its properties relative to the AL distribution. In Section 3, we formulate the Bayesian quantile regression model, including a prior spec- 3 ification for the regression coefficients that encourages shrinkage resulting in regularized quantile regression, andaTobit quantileregression formulation. InSection 4, we presentresultsfrom asim- ulation study to compare the performance of the AL and the proposed distribution in regularized quantile regression. The methodology is illustrated with three data examples in Section 5, focusing again on comparison with the AL quantile regression model. Finally, Section 6 concludes with a summaryanddiscussionofpossibleextensions. 2 The generalized asymmetric Laplace distribution Theconstructionofthenewdistributionismotivatedbythemostcommonlyusedmixturerepresen- tationoftheALdensity. Inparticular, (cid:90) fAL(y | µ,σ) = N(y | µ+σA(p)z,σ2B(p)z)Exp(z | 1)dz (2) p R+ whereA(p) = (1−2p)/{p(1−p)}andB(p) = 2/{p(1−p)}. Moreover,N(m,W)denotesthenor- maldistributionwithmeanmandvarianceW,andExp(1)denotestheexponentialdistributionwith mean1. Weusesuchnotationthroughouttoindicateeitherthedistributionoritsdensity,depending onthecontext. Themixtureformulationin(2)enablesexplorationofextensionstotheALdistribution. Extend- ing the Exp(1) mixing distribution is not a fruitful direction in terms of evaluation of the intergal, and, moreimportantly, withrespecttofixingpercentilesoftheresultingdistribution. However, both goals are accomplished by replacing the normal kernel in (2) with a skew normal kernel (Azzalini, 1985). In its original parameterization, the skew normal density is given by fSN(y | ξ,ω,λ) = 2ω−1φ(ω−1(y−ξ))Φ(λω−1(y−ξ)),whereφ(·)andΦ(·)denotethedensityanddistributionfunc- tion, respectively, of the standard normal distribution. Here, ξ ∈ R is a location parameter, ω > 0 a scaleparameter,andλ ∈ Rtheskewnessparameter. Keytoourconstructionisthefactthattheskew normaldensitycanbewrittenasalocationnormalmixturewithmixingdistributiongivenbyastan- dardnormaltruncatedonR+ (Henze,1986). Morespecifically,reparameterize(ξ,ω,λ)to(ξ,τ,ψ), where τ > 0 and ψ ∈ R, such that λ = ψ/τ and ω = (τ2 + ψ2)1/2. Then, fSN(y | ξ,τ,ψ) = 4 (cid:82) N(y | ξ + ψs,τ2)N+(s | 0,1)ds, where N+(0,1) denotes the standard normal distribution R+ truncatedoverR+. The proposed model, referred to as generalized asymmetric Laplace (GAL) distribution, is built byaddingashapeparameter,α ∈ R,tothemeanofthenormalkernelin(2)andmixingwithrespect to a N+(0,1) variable. More specifically, the full mixture representation for the density function, f(y | p,α,µ,σ),ofthenewdistributionisasfollows (cid:90)(cid:90) N(y | µ+σαs+σA(p)z,σ2B(p)z)Exp(z | 1)N+(s | 0,1)dzds. (3) R+×R+ Note that, integrating over s in (3), the GAL density can be expressed in the form of (2) with the N(y | µ+σA(p)z,σ2B(p)z) kernel replaced with a skew normal kernel, which, in its original parameterization,haslocationparameterµ+σA(p)z,scaleparameterσ{α2+B(p)z}1/2,andskewness parameterα{B(p)z}−1/2. Evidently,whenα = 0,f(y | p,0,µ,σ)reducestotheALdensity. To obtain the GAL density, we integrate out first z and then s in (3). The integrand of (cid:82) N(y | µ + σαs + σA(p)z,σ2B(p)z)Exp(z | 1)dz can be recognized as the kernel of a R+ generalized inverse-Gaussian density. Therefore, integrating out z, we obtain f(y | p,α,µ,σ) = (cid:82) p(1−p)σ−1exp(cid:8)−σ−1[p−I(y < µ+σαs)][y−(µ+σαs)](cid:9) N+(s | 0,1)ds. Thisintegral R+ involves a normal density kernel, but care is needed with the limits of integration which depend on the sign of y −µ and of α. Combining the resulting expressions from all possible cases, we obtain thatforα (cid:54)= 0,theGALdensityisgivenby p(1−p) (cid:18)(cid:20) (cid:18)y∗ (cid:19) (cid:21) (cid:26) 1 (cid:27) (cid:18)y∗ (cid:19) f(y |p,α,µ,σ) = 2 Φ −p α −Φ(−p α) exp −p y∗+ (p α)2 I >0 σ α α− α− α− 2 α− α (cid:20) y∗ (cid:18)y∗ (cid:19)(cid:21) (cid:26) 1 (cid:27)(cid:19) +Φ p α− I >0 exp −p y∗+ (p α)2 (4) α+ α α α+ 2 α+ wherey∗ = (y−µ)/σ,p = p−I(α > 0),p = p−I(α < 0),withp ∈ (0,1). Therelatively α+ α− complexformofthedensityin(4)isnotanobstaclefromapracticalperspective,sinceitshierarchical mixturerepresentationfacilitatesstudyofmodelpropertiesandMarkovchainMonteCarloposterior simulation. There is a direct link between the GAL distribution and the p th quantile for any p ∈ (0,1); 0 0 5 note that parameter p no longer corresponds to the cumulative probability at the quantile for α (cid:54)= 0. (cid:82)µ When α > 0, the distribution function of (4) at µ is given by f(y | p,α,µ,σ)dy = 2pΦ[(p− −∞ 1)α]exp(cid:8)(p−1)2α2/2(cid:9). Hence,lettingγ = (1−p)α,thedistributionfunctionbecomes, (cid:90) µ f(y | p,γ,µ,σ)dy = pg(γ) with g(γ) = 2Φ(−|γ|)exp(γ2/2). −∞ Weuse|γ|above,sincethisisthegeneralformofg(γ)thatappliesalsointheα < 0case. Note that, for γ ∈ R−, dg(γ)/dγ = 2h(γ)exp(γ2/2), where h(γ) = φ(γ) + γΦ(γ). The function h(γ) is monotonically increasing in R−, since dh(γ)/dγ = Φ(γ) > 0. Moreover, h(0) = (2π)−1/2 > 0, and lim h(γ) = 0. Therefore, h(γ) > 0 for γ ∈ R−, and thus g(γ) is γ→−∞ monotonicallyincreasinginR−. Sinceg(γ)isanevenfunction,italsoobtainsthatitismonotonically decreasinginR+. (cid:82)µ Considernowsetting f(y | p,γ,µ,σ)dy =pg(γ) = p . Then,thefactthatg(γ)isdecreas- −∞ 0 inginR+ combinedwithg(γ) > p , implythatforeachγ > 0inthedomainthatrespectsthecon- 0 (cid:82)µ ditionofp ∈ (0,1)andα > 0,thereisauniquesolutionofpthatensures f(y | p,γ,µ,σ)dy = −∞ (cid:82)∞ p ,andsubsequentlyauniqueαbasedonγ = (1−p)α. Forα < 0,setting f(y | p,γ,µ,σ)dy = 0 µ 1−p andlettingγ =pαleadstothesameargument. 0 Theaboveconnectionbetween(p ,γ)and(p,α)suggeststhatbyreparameterizationwithdesired 0 p andγ = [I(α > 0)−p]|α|,wecanderiveanewfamilyofdistributionswiththepercentileforfixed 0 p givenbyµ,andwithanadditionalshapeparameterγ. Forγ (cid:54)= 0,thedensity,f (y | γ,µ,σ),of 0 p0 suchquantile-fixedGALdistributionis p(1−p)(cid:32)(cid:26) (cid:18) p y∗ p (cid:19) (cid:18)p (cid:19)(cid:27) (cid:40) γ2 (cid:18)p (cid:19)2(cid:41) (cid:18)y∗ (cid:19) 2 Φ − γ+ + γ−|γ| −Φ γ−|γ| exp −p y∗+ γ− I >0 σ |γ| p p γ− 2 p γ γ+ γ+ γ+ (cid:20) p y∗ (cid:18)y∗ (cid:19)(cid:21) (cid:26) γ2(cid:27)(cid:19) +Φ −|γ|+ γ+ I >0 exp −p y∗+ (5) |γ| γ γ+ 2 where p ≡ p(γ,p ) = I(γ < 0) + {[p − I(γ < 0)]/g(γ)}, p = p − I(γ > 0), p = 0 0 γ+ γ− p−I(γ < 0),andy∗ =(y−µ)/σ. Parameterγ hasboundedsupportoverinterval(L,U),whereL is the negative root of g(γ) = 1−p and U is the positive root of g(γ) = p . For instance, γ takes 0 0 valuesin(−0.07,15.90),(−1.09,1.09)and(−2.90,0.39)whenp = 0.05,p = 0.5andp = 0.75, 0 0 0 6 p =0.05 p =0.5 p =0.75 0 0 0 0 5 2 8 g= 0 0.2 g= 0 0. g= 0 0.0 g= −0.04 g= −0.8 g= −1.5 g= 4 20 g= −0.4 5 g= −0.5 6 g= 8 0. g= 0.6 0.1 g= 0.2 0 0. 5 Density 0.04 Density 100.1 Density 0.10 0. 5 0.02 0.05 0.0 0 0 0 0 0 0 0. 0. 0. −20 −10 0 10 20 −10 −5 0 5 10 −10 −5 0 5 10 e e e Figure1: Densityfunctionofquantile-fixedgeneralizedasymmetricLaplacedistributionwithµ = 0, σ = 1 anddifferentvaluesofγ,forp =0.05,0.5and0.75. Inallcases,thesolidlinecorrespondstotheasymmetric 0 Laplacedensity(γ =0). respectively. Whenγ = 0,thedensityreducestotheALdensity,whichisalsoalimitingcaseof(5). Thedensityfunctioniscontinuousforallpossibleγ values. The quantile-fixed GAL distribution has three parameters, µ, σ and γ. Note that Y has density f (· | γ,µ,σ) if and only if (Y − µ)/σ has density f (· | γ,0,1). Hence, similarly to the AL p0 p0 distribution,µisalocationparameterandσisascaleparameter. Thenewshapeparameterγ enables the extension relative to the quantile-fixed AL distribution. As demonstrated in Figure 1, γ controls skewness and tail behaviour, allowing for both left and right skewness when the median is fixed, as wellasforbothheavierandlightertailsthantheasymmetricLaplace,thedifferencebeingparticularly emphatic for extreme percentiles. Moreover, as γ varies, the mode is no longer held fixed at µ; it is less than µ when γ < 0 and greater than µ when γ > 0. The above attributes render the proposed distributionsubstantiallymoreflexiblethantheALdistribution. Finally, we note that parameter γ satisfies likelihood identifiability. Consider the location-scale standardized density, f (· | γ,0,1), which is effectively the model for the errors in quantile regres- p0 sion. Then, assume f (y | γ ,0,1) = f (y | γ ,0,1), for all y ∈ R. Given that parameter γ p0 1 p0 2 7 controlsthemodeofthedensity, thisimpliesthatγ andγ musthavethesamesign. Workingwith 1 2 eitherofthetwocases(thatis,γ > 0andγ > 0orγ < 0andγ < 0)inexpression(5),wearrive 1 2 1 2 atg(γ ) =g(γ ),which,basedonthemonotonicityoffunctiong(·),impliesγ = γ . 1 2 1 2 3 Bayesian quantile regression with GAL errors 3.1 Inferenceforlinearquantileregression Considercontinuousresponsesy andtheassociatedcovariatevectorsx ,fori = 1,...,n. Thelinear i i quantileregressionmodelissetupasy =xTβ+(cid:15) ,wherethe(cid:15) ariseindependentlyfromaquantile- i i i i (cid:82)0 fixedGALdistributionwith f ((cid:15) | γ,0,σ)d(cid:15) =p . Owingtothemixturerepresentationofthe −∞ p0 0 newdistribution,themodelforthedatacanbeexpressedhierarchicallyasfollows y | β,γ,σ,z ,s in∼d. N(y | xTβ+σC|γ|s +σAz ,σ2Bz ), i = 1,...,n i i i i i i i i z ,s in∼d. Exp(z | 1)N+(s | 0,1), i = 1,...,n (6) i i i i whereC = [I(γ > 0)−p]−1,andAandB arethefunctionsofpgivenin(2). Sincepisafunction of γ and p , A, B and C are all functions of parameter γ. The Bayesian model is completed with 0 priorsforβ,σ andγ. Here,weassumeanormalpriorN(m ,Σ )forβ andaninverse-gammaprior 0 0 IG(a ,b )forσ,withmeanb /(a −1)provideda > 1. Foranyspecifiedp ,γ isdefinedoveran σ σ σ σ σ 0 interval (L,U) with fixed finite endpoints, and thus a natural prior for γ is given by a rescaled Beta distribution,withtheuniformdistributionavailableasadefaultchoice. The augmented posterior distribution, which includes the z and the s , can be explored via a i i MarkovchainMonteCarloalgorithmbasedonGibbssamplingupdatesforallparametersotherthan γ. AsinKozumiandKobayashi(2011),wesetv =σz ,i = 1,...,n. Then,theposteriorsimulation i i methodisbasedonthefollowingupdates. 1. Sampleβ fromN(m∗,Σ∗),withcovariancematrixΣ∗ =[Σ−1+(cid:80)n x xT/(Bσv )]−1 and 0 i=1 i i i meanvectorm∗ =Σ∗{Σ−1m +(cid:80)n x [y −(σC|γ|s +Av )]/(Bσv )}. 0 0 i=1 i i i i i 2. Foreachi = 1,...,n,samplev fromageneralizedinverse-Gaussiandistribution,GIG(0.5,a ,b ), i i i 8 wherea = [y −(xTβ+σC|γ|s )]2/(Bσ)andb = 2/σ+A2/(Bσ),withdensitygivenby i i i i i GIG(x | ν,a,b) ∝xν−1exp{−0.5(a/x+bx)}. 3. Foreachi = 1,...,n,samples fromanormalN(µ ,σ2)distributiontruncatedonR+,where i si si σ2 =[(Cγ)2σ/(Bv )+1]−1 andµ =σ2C|γ|[y −(xTβ+Av )]/(Bv ). si i si si i i i i 4. Sampleσ fromaGIG(ν,c,d)distribution, whereν =−(a +1.5n), c =2b +2(cid:80)n v + σ σ i=1 i (cid:80)n [y −(xTβ+Av )]2/(Bv ),andd = (cid:80)n (Cγs )2/(Bv ). i=1 i i i i i=1 i i 5. UpdateγwithaMetropolis-Hastingstep,usinganormalproposaldistributiononthelogitscale over(L,U). Basedonthehierarchicalmodelstructure,theposteriorpredictiveerrordensitycanbeexpressed as p((cid:15) | data) = (cid:82) N((cid:15) | σC|γ|s+σAz,σ2Bz)Exp(z | 1)N+(s | 0,1)π(γ,σ|data)dsdzdγdσ, andthusestimatedthroughMonteCarlointegration,usingtheposteriorsamplesof(γ,σ). 3.2 Quantileregressionwithregularization Since the GAL distribution is constructed through modifying the mixture representation of the AL distribution,itretainssomeoftheinterestingpropertiesoftheALdistribution. Inparticular,working with the hierarchical representation of the GAL distribution, we are able to retrieve an extended versionofthechecklossfunctionwhichcorrespondstoasymmetricLaplaceerrors. Consider the collapsed posterior distribution, π(β,γ,σ,s ,...,s | data), that arises from (6) by 1 n marginalizingoverthez . Then, thecorrespondingposteriorfullconditionalforβ canbeexpressed i as (cid:40) n (cid:41) 1 (cid:88) π(β | γ,σ,s ,...,s ,data) ∝ π(β)exp − ρ (y −xTβ−σH(γ)s ) 1 n σ p i i i i=1 where π(β) is the prior density for β, H(γ) = γg(γ)/{g(γ)−|p −I(γ < 0)|}, and p = I(γ < 0 0)+{[p −I(γ < 0)]/g(γ)},withp theprobabilityassociatedwiththespecifiedquantilemodeled 0 0 through xTβ. Hence, ignoring the prior contribution, finding the mode of the posterior full condi- i tional for β is equivalent to minimizing with respect to β the adjusted loss function (cid:80)n ρ (y − i=1 p i 9 xTβ −σH(γ)s ); note that in the special case with asymmetric Laplace errors, that is, for γ = 0, i i thisreducestothechecklossfunctionwithp =p . 0 Based on the above structure, the positive-valued latent variables s can be viewed as response- i specific weights that are adjusted by real-valued coefficient H(γ), which is fully specified through theshapeparameterγ. Theresultisthereal-valued,response-specifictermsσH(γ)s ,whichreflect i ontheestimationofβtheeffectofoutlyingobservationsrelativetotheALdistribution. Apromising directiontofurtherexplorethisstructureisinthecontextofvariableselection. Forinstance,Lietal. (2010) study connections between different versions of regularized quantile regression and different priorsforβ, workingwithasymmetricLaplaceerrors. Themainexampleislassoregularizedquan- tile regression, which can be connected to the Bayesian asymmetric Laplace error model through a hierarchical Laplace prior for β. We consider this prior below extending the AL error distribution to the proposed GAL distribution. The perspective we offer may be useful, since it can be used to exploreregularizationadjustingthelossfunction,throughtheresponsedistribution,inadditiontothe penaltyterm,throughthepriorfortheregressioncoefficients. Here,wedenotebyβ thed-dimensionalvectorofregressioncoefficientsexcludingtheintercept β . Then,theLaplaceconditionalpriorstructureforβ isgivenby 0 (cid:89)d λ (cid:26) λ (cid:27) (cid:89)d (cid:90) 1 (cid:26) β2 (cid:27) η2 (cid:26) η2 (cid:27) π(β | σ,λ) = exp − |β | = √ exp − k exp − ω dω . k k k 2σ σ R+ 2πωk 2ωk 2 2 k=1 k=1 ThesecondexpressionaboveutilizesthenormalscalemixturerepresentationfortheLaplacedistribu- tion,whichhasbeenexploitedforposteriorsimulationinthecontextoflassomeanregression(Park and Casella, 2008). Moreover, to facilitate Markov chain Monte Carlo sampling, we reparameterize in terms of η = λ/σ and place a gamma prior on η2. The lasso regularized version of model (6) is completedwithanormalpriorforβ ,andwiththepriorsfortheotherparametersasgiveninSection 0 3.1. The posterior simulation algorithm is the same with the one described in Section 3.1 with the exceptionoftheupdatesfortheβ ,k = 1,...,d,andforη2. Usingthemixturerepresentationofthe k Laplaceprior,eachβ canbesampledfromanormaldistribution,whereasη2 hasagammaposterior k fullconditionaldistribution. 10

A New Family of Error Distributions for Bayesian Quantile Regression

Similar A New Family of Error Distributions for Bayesian Quantile Regression