Bayesian Estimation of Simultaneous Regression Quantiles Using Hamiltonian Monte Carlo

: The simultaneous estimation of multiple quantiles is a crucial statistical task that enables a thorough understanding of data distribution for robust analysis and decision-making. In this study, we adopt a Bayesian approach to tackle this critical task, employing the asymmetric Laplace distribution (ALD) as a flexible framework for quantile modeling. Our methodology implementation involves the Hamiltonian Monte Carlo (HMC) algorithm, building on the foundation laid in prior work, where the error term is assumed to follow an ALD. Capitalizing on the interplay between two distinct quantiles of this distribution, we endorse a straightforward and fully Bayesian method that adheres to the non-crossing property of quantiles. Illustrated through simulated scenarios, we showcase the effectiveness of our approach in quantile estimation, enhancing precision via the HMC algorithm. The proposed method proves versatile, finding application in finance, environmental science, healthcare, and manufacturing, and contributing to sustainable development goals by fostering innovation and enhancing decision-making in diverse fields.


Introduction
Quantile regression, which is garnering increased attention from both theoretical and empirical perspectives [1][2][3][4][5], serves as a statistical cornerstone for modeling the intricate relationship between a set of predictor variables and a response variable [6,7].This method excels in estimating the conditional quantiles of the response variable, going beyond conventional mean estimation.By shedding light on the varying impacts of predictor variables across different segments of a distribution, quantile regression offers a nuanced understanding of complex relationships in data analysis.As researchers delve into unraveling the intricacies of variable interplay, the application of quantile regression emerges as an invaluable approach for capturing and interpreting diverse patterns.For instance, we can enumerate examples in various fields such as the economy [8,9], ecology [10], epidemiology [11], survival analysis [12], and environmental economics [13,14].Furthermore, quantile regression is a form of statistical and econometric regression analysis that adds a particular method to estimating families of conditional quantile functions and supplements the exclusive focus of least squares techniques based on the estimation of conditional mean functions [15].Originating from the seminal work of [16], quantile regression introduces a versatile method for estimating conditional quantiles.
The authors demonstrated that any linear quantile of order τ could be derived from minimizing a check function, ρ τ (u) = u(τ − 1 u<0 ), where 1 represents the indicator function.In exclusively linear regression models, Koenker et al. [16] suggest a straightforward method to estimating conditional quantiles, consider the regression linear normal model, as the solution of the minimization problem that follows: where x t , t = 1, . . ., T denote a sequence of (row), x ′ is the transposed of x, K-vector of a known design matrix, and y t , t = 1, . . ., T are a random sample on the regression process u = y − x ′ t β having distribution function F.Then, the τ-th sample quantile, 0 < τ < 1, may be defined as any solution of the minimization problem (1).
In the framework of statistical research, Bayesian approaches have emerged as effective tools for investigating quantile regression [17][18][19][20] after the notable developments presented in [21].The authors utilize a likelihood function grounded in the asymmetric distribution of Laplace to advance the understanding of quantile regression.While both frequentist and Bayesian perspectives prove suitable for examining individual quantiles [22], the complexities of real-world applications often demand more nuanced exploration.To intricately capture the relationship between the response variable and explanatory variables, it becomes imperative to estimate regression quantiles at various orders.Although it is conceivable to perform separate estimations for these quantiles, this approach proves unsatisfactory; particularly in scenarios involving multiple explanatory variables, it provides estimators that do not respect the property of the non-crossing of the regression quantiles.
To circumvent this limitation, a more favorable approach involves the simultaneous estimation of multiple regression quantiles [23][24][25][26].This ensures a coherent and consistent estimation procedure, particularly in situations characterized by the presence of multiple explanatory variables.
In pursuit of refining the Bayesian quantile estimation paradigm, this article directs its attention to the groundwork laid by Merhi Bleik in the seminal work presented in [27].The author's utilization of the Metropolis-Hastings within Gibbs algorithm proved instrumental in sampling unknown parameters from their full conditional distribution, offering a new methodological framework.The proposed estimation procedure, SBQR, for simultaneous Bayesian quantile regression, as presented in [27], guarantees the fundamental property of non-crossing, under the assumption that the Asymmetric Laplace distribution (ALD) is the underlying data distribution.In extending and enhancing this foundational approach, our primary objective is to present a unified computational methodology.This methodology is designed to facilitate the simultaneous estimation of regression quantiles, employing the advanced Hamiltonian Monte Carlo (HMC) algorithm.Building on the achievements of [27], our proposed methodology aims to overcome potential limitations and further streamline the quantile estimation process.By integrating the HMC algorithm, we seek to harness its efficiency in exploring parameter spaces, thus enhancing the precision of quantile estimates.To evaluate the performance of our algorithm, we compare the mean of the Bayesian estimators with the theoretical quantile values using the root mean square error.For this reason, our objective is to simultaneously estimate the Bayesian regression quantile using the HMC algorithm.The main contributions of the paper are then the following:

•
Conduct a simultaneous quantile estimation embedded with the HMC algorithm through Bayesian inference; • Conduct a sensitivity analysis on the performance of the HMC algorithm in handling univariate linear quantiles and multivariate linear quantiles.
This article serves as a comprehensive guide to the unified computational framework, providing researchers and practitioners with a robust tool for simultaneous regression quantile estimation.By addressing challenges associated with multiple quantile modeling in a cohesive manner, our methodology proves versatile across diverse fields when used with ALD data.Notably, applications in finance, where risk assessment demands a nuanced understanding of the distribution, benefit from the precision offered by simultaneous regression quantile estimation.In environmental science, the ability to capture extreme events accurately is crucial, while in healthcare, predicting patient outcomes relies on a comprehensive grasp of varying response scenarios.Additionally, in manufacturing, ensuring quality control is facilitated by an in-depth understanding of distributional patterns.This versatility aligns with the broader societal goals, contributing to Sustainable Development Goal 9: Industry, Innovation, and Infrastructure.By fostering innovation in statistical methodologies and decision-making processes, our approach seeks to enhance the well-being of societies and align with broader sustainability objectives.
This paper introduces a novel approach by combining existing research with the innovative application of the HMC algorithm to simultaneous quantile estimation, offering significant advancements in the field.The paper is organized as follows: The ALD and its Gaussian mixture are presented in Section 2, in addition to the method proposed by the author in [27] for simultaneous Bayesian quantiles.A brief overview of HMC is presented in Section 3. The performance of the HMC algorithm is examined in cases of univariate linear quantiles and multivariate linear quantiles on a simulated data set in Section 4. Finally, Section 5 highlights the paper's contributions and provides a discussion of perspectives.

Asymmetric Laplace Distribution
The ALD is a probability distribution that is used to model data with asymmetric tails and skewness.It is a continuous probability distribution defined on the real line.The distribution is characterized by three parameters: location, scale, and skewness parameters.Indeed, the ALD is often used in Bayesian statistics, particularly in the context of Bayesian regression models that involve quantile regression.
A random variable U follows an ALD(µ, σ, τ) if its probability density is given by where 0 < τ < 1 is the skew parameter, σ > 0 is the scale parameter, and −∞ < µ < ∞ is the location parameter.The mean and variance of the ALD with density (2) are Some other ALD properties can be found in [28].In order to simplify the posterior inference of the quantiles, it is convenient to use the location scale mixture representation.In particular, we are obliged to use it in HMC because ALD is not differentiable and to calculate the full conditional distributions (for more details, see [20]).The location-scale mixture of the ALD distribution will be used.Let ω be an exponential latent variable with parameter 1/σ denoted by exp(1/σ) and let Z be a standard normal variable so that ω and Z are independent, then ϵ ∼ ALD(0, σ, τ) has the following representation: where

Regression Quantile
We begin with a fundamental definition of sample quantiles, a concept readily extended to the linear regression model without the typical need for an ordered set of sample observations.Consider the set y t : t = 1, . . ., T of independent and identically distributed (i.i.d.) observations on a random variable Y governed by the distribution function F. The τ-th sample quantile, where 0 < τ < 1, can then be characterized as any solution to the minimization problem: where ρ τ is the check loss function, defined as (3) Equivalently, we may write the minimization problem as min β∈R [ ∑ t∈{t:y t ⩾β} Huber et al. [29] observed the difficulty in identifying outliers within the context of regression, emphasizing the inherent confusion in extending standard notions of sample quantiles to linear models based on the ordering of sample observations.This confusion is effectively addressed by a direct generalization of the earlier-described minimization problem.Consider a sequence of k-vectors from a known design matrix denoted by {x t , t = 1, . . ., T}. Assume that {y t , t = 1, . . ., T} constitutes a random sample from the regression process, where u t = y t − x ′ t β follows a distribution function F such that F(0) = τ.The τ-th quantile of the linear regression model, with 0 < τ < 1, is defined as any solution to the minimization problem:

Bayesian Linear Quantile Regression
Consider the following quantile regression model: We assume that the error ϵ and x are independent, and ϵ ∼ F, F is unknown, but F(0) = τ.Hence, x ′ β τ denotes the τ-th quantile of Y|x.For all x i ∈ R, suppose the relationship between ι τ (x i ) and x i could be modeled as ι τ (x i ) = x i β τ , where β τ is a vector with unknown interest parameters.From i.i.d (x t , Y t ) t=1,...,T issued from the model defined in Equation ( 5), we are interested in estimating quantiles of the Y|x law via the Bayesian approach.It then suffices to estimate β τ .If F is unknown, the likelihood is not available, but the Bayesian approach requires some likelihood; hence, we use the "working likelihood" defined as the ALD: where µ = x ′ t β τ and ρ τ (u) is as defined in (3).The τ-th quantile of this distribution is zero; hence, the pseudo posterior is Equation ( 6) describes the ALD probability distribution, where y t is the observed data point, x ′ t is the design matrix for observation t, β is a set of parameters, σ is a scale parameter, and τ is a quantile parameter.This equation provides a natural way to solve the problem of Bayesian estimation of the quantile q τ of order τ.Moreover, the relationship between minimizing the loss function (3) and maximizing a likelihood function is demonstrable.This equivalence arises when the likelihood function is constructed by combining independently distributed asymmetric densities of the Laplace distribution.Furthermore, it is noteworthy to emphasize that the simple Bayesian quantile can be regarded as a specific instance of a Bayesian regression quantile.In this particular case, the simplicity is characterized by setting k = 1 and considering only the variable x, with the absence of additional co-variables.This distinction highlights a specialized scenario within the broader framework of Bayesian quantile regression.

Simultaneous Quantile Regression
Many approaches have been proposed to overcome the crossing problem.See, among others, in the frequentist context [30][31][32].In the simultaneous Bayesian context, on the other hand, Ref. [33] suggested a two-stage semi-parametric regression method.Steven G Xu and Brian J Reich [34] have also proposed a novel treatment to nonlinear simultaneous quantile regression.They have achieved this by specifying a Bayesian nonparametric model for the conditional distribution.In a similar vein, Das and Ghosal [35] represent the quantile process by articulating it as a weighted sum of B-spline basis functions corresponding to the quantile level.In this study, the method proposed in [27] is considered, which assumes the error term has an ALD.It is noteworthy that imposing a given distribution on u t is restrictive but provides a fully Bayesian approach to handling simultaneous regression quantiles of different orders.This is performed first by partitioning the entire sample into s sub-samples, where s is the number of quantiles of interest; second, by using the relationship between any two ALD quantiles; and third, by rewriting the entire likelihood into s terms where the j-th term depends only on q τ j , where q τ j is the j-th term of the vector of quantiles.Let us explain how it works: 1.
Consider partitioning the entire sample into the sub-samples below: where I j = {(j − 1)r + 1, . . ., jr} and with assume that r = n s is an integer without loss of generality.

2.
In order to characterize the likelihood through all quantiles of interest, we use the relationship that connects any τ-th quantile to the p-th quantile of the ALD(q p , σ, p): where

3.
We rewrite the model given by (5) from Equation ( 7) as follows: where q τ j ,I j = (q τ j (x i j )) i j ∈I j , g(τ j , p) = g(τ j , p)1 i j ∈I j i j ∈I j , The likelihood connected with the model provided by Equation ( 5) is then rewritten as the product of s likelihoods, each depending on only one sub-sample: L(q τ 1 ,I 1 , . . . ,q τ s ,I s , σ, p; x I 1 , y I 1 , . . . ,

Hamiltonian Monte Carlo
Markov chain Monte Carlo (MCMC) was employed to simulate state distributions for an idealized molecular system.An alternative molecular dynamics approach was implemented, featuring deterministic molecular movement adhering to Newton's law of motion elegantly formalized as Hamiltonian dynamics.The inaugural MCMC algorithm incorporating Hamiltonian dynamics was the static HMC [36].HMC [37,38] has favorable properties in complex models.The sampling from Metropolis-Hastings [39] and Gibbs [40] relies on random samples from an easy-to-sample distribution of proposals or conditional densities.In complex problems, these algorithms may not give good results; for example, random walks have difficulty "exploring" distributions.The purpose of the HMC algorithm is to use a different kind of sampling that is expected to more effectively "explore" the posterior distribution.Moreover, the HMC requires the differentiability of the density function.According to [41], HMC's strength is its ability to explore high-dimensional parameter spaces more effectively.By utilizing gradient information from the posterior distribution, HMC avoids the random-walk nature of traditional MCMC methods.This results in a faster convergence, lower autocorrelation, and enhanced mixing, which provides more accurate estimates of the posterior distribution using fewer samples (a simple example is presented in Appendix B).In addition, HMC outperforms Gibbs sampling and Metropolis-Hastings by efficiently exploring high-dimensional target distributions without succumbing to random walk behavior or sensitivity to correlated parameters [42].Beyond its ability to manage complex and multimodal distributions, HMC's versatility extends to diverse fields, making it a useful tool in areas such as machine learning, physics, and more.
HMC, essentially a variant of Metropolis-Hastings, draws inspiration from a physical analogy with a distinct approach to proposal generation.The dynamics can be visualized in two dimensions, such as that of a frictionless puck that slides across a variable-height surface.The state of this system comprises the puck position, provided by a two-dimensional vector z, and the puck momentum, provided by a two-dimensional vector p (its mass times its speed).The puck's potential energy, U(z), is proportional to its present position's ground height, while its kinetic energy, K(p), is equal to p 2 /(2m), where m is the puck's mass.The puck passes at a steady velocity, equivalent to p/m, at a level part of the surface.The potential energy for these factors forms the negative logarithm of the probability density.
This physical analogy serves as a conceptual tool for understanding the dynamic interplay between position and momentum in the system, providing valuable insight into the foundational principles of the HMC algorithm.Subsequently, in HMC applications, Hamiltonian functions are commonly employed and can be represented as follows: Here, U(z) is called potential energy and is defined as minus the distribution log probability density for z (U(z) = − log(l(x|z) × g(z)) that we want to sample, plus any useful constant.K(p) is called kinetic energy and is generally described as Here, M is a symmetric, positive-definite "mass matrix".After identifying our potential and kinetic energy functions and introducing momentum variables that are relevant to our position variables, we can simulate a Markov chain in which each iteration takes place:

•
Randomly sample K(p) momentum variables, regardless of present position; • Update Metropolis using simulated Hamiltonian dynamics (implemented by the leapfrog method, for example).
More specifically, the momentum variable determines where the position z passes during this second dynamic proposal, and the gradient of U(z) determines how the momentum p changes in accordance with the equations of Hamilton: Some properties of Hamiltonian dynamics are reversibility, conservation of the Hamiltonian, volume preservation, and symplecticness (for more details, see [38]).
Leapfrog Method: Hamilton's equations need to be approximated by discretizing time for software execution using some tiny step size, ε.We do this by assuming that the Hamiltonian has the form H(z, p) = U(z) + K(p), as in Equation ( 8), and by assuming that kinetic energy has the form K(p) = p T M −1 p/2, as in Equation ( 9); moreover, M is diagonal, with m 1 , . . ., m d diagonal components, so that The leapfrog method works as follows: 1. 2.
For the momentum variables, we begin with a half step, then we perform a complete step for the position variables using the current momentum variables' values, and lastly, we perform another half step for the momentum variables using the current position variables' values.

HMC Algorithm:
Each iteration of the HMC algorithm has two steps: 1. New values are drawn randomly from their Gaussian distribution for the momentum variables, irrespective of the current position variables values.

2.
Updating Metropolis is performed using Hamiltonian dynamics to suggest a new state.Starting with the current state (z, p), Hamiltonian dynamics is simulated using the leapfrog method (or some other reversible method that preserves volume) for L steps, using a step size of ε.At the end of this L-step trajectory, the momentum variables are then negated, resulting in a suggested state (z * , p * ).This proposed state is accepted as the next state in the Markov chain with probability: If the proposed state is not approved (i.e., dismissed), the following state will be the same as the present state.

Numerical Results
In this section, we will perform a Bayesian quantile regression procedure to estimate simultaneous regression quantiles using HMC.HMC will be used through the probabilistic programming language STAN [43][44][45] program to estimate simultaneous linear quantiles.STAN provides a flexible framework for specifying and fitting a wide range of Bayesian models.Numerous probability densities, matrix operations, and numerical ODE solvers are supported by the expressive language STAN [46].
Two different models will be considered for the τ-th quantile for the model provided by ( 5): 1.
Our methodology involves the utilization of the state-of-the-art STAN program, leveraging Hamiltonian Monte Carlo for statistical modeling and high-performance statistical computing.This platform has been extensively employed across diverse disciplines, including the social, biological, and physical sciences, engineering, and industry, for tasks such as statistical modeling, data analysis, and prediction.Given the inherent relationship between any two ALD quantiles, we recognize that distinct quantiles vary solely in terms of their intercepts β τ j .Consequently, all quantiles of interest have the same slope ; then, the unknown vector parameter is β.The prior of β is Gaussian, with mean 0 and variance 1.For further details about the STAN code, refer to the Appendix A.

Univariate Linear Quantile
Suppose the error terms ϵ i=1,...,T i.i.d ∼ ALD(0, 0.1, 0.75).As the Bayesian approach requires a likelihood, we use the ALD to make the Bayesian inference, and our quantiles of interest are τ = 0.25 and τ = 0.5; hence, we use the relationship between the two ALDs mentioned in Section 2.
To study the quality of the estimate of the quantiles, we calculate the root mean square error: where q τ is the theoretical quantile and qτ is the posterior mean quantile regression.
To assess the impact of data variance on estimation, we assumed that the error terms ϵ i=1,...,T i.i.d ∼ ALD(0, 0.3, 0.75).The obtained results are presented in Table 1.We run our algorithm for 20,000 iterations, with half of them being burn-in.Gelman-Rubin diagnostic (R-hat) values are depicted in Figure 1, with R approaching 1, indicating the convergence of multiple chains.The densities of the estimated βs are depicted in Figure 2, providing insight into the distribution patterns.
Furthermore, Figure 3 illustrates that the autocorrelation of each chain converges to 0. While comparing the HMC method used in this paper with the Metropolis-Hastings within Gibbs algorithm used in [27], we find that using HMC exhibits lower autocorrelation at earlier lags (around 0 after 10 lags) compared to the Metropolis-Hastings within Gibbs algorithm, where autocorrelation remains around 0 only after 50 lags.This suggests that HMC achieves greater independence between samples sooner, indicating potentially more efficient sampling and better exploration of the distribution space compared to the algorithm used in [27].
Additionally, the examination of the trace plots in Figure 4 reveals consistent and stable behavior in the MCMC chains, showcasing no discernible oscillation, drift, or other abnormalities.This pattern suggests a thorough exploration of the parameter space, indicating successful convergence to the target distribution.The small variance observed in Table 1 further supports the reliability of the results.Furthermore, the Monte Carlo standard error (MCSE) is 0 for all βs, indicating that the estimated values are likely close to the true values, as shown in Table 1.In addition, when the variance of the data is larger, the estimation worsens (higher RMSE and higher variance).This is normal, as a larger variance makes it harder to estimate the parameters.
It should be noted that diagnostics can only reliably be used to determine a lack of convergence and not detect convergence [47].Therefore, from the obtained results, the trace plots do not seem to indicate non-stationary behavior.

Multivariate Linear Quantile
Suppose the error terms ϵ i=1,...,T i.i.d ∼ ALD(0, 0.1, 0.75).As the Bayesian approach requires a likelihood, we use the ALD to make the Bayesian inference, and the quantiles of interest are τ = 0.1 and τ = 0.12; hence, we use the relationship between the two ALDs mentioned in Section 2.
The Gelman-Rubin diagnostic (R-hat) values are depicted in Figure 5, with R approaching 1, indicating the convergence of multiple chains.The histograms of the estimated βs are depicted in Figure 6, providing insights into the distribution patterns.
Furthermore, Figure 7 illustrates that the autocorrelation of each chain converges to 0. While comparing the HMC method used in this paper with the Metropolis-Hastings within Gibbs algorithm used in [27], we find that using HMC exhibits lower autocorrelation at earlier lags (around 0 before 10 lags) compared to the Metropolis-Hastings within Gibbs algorithm, where autocorrelation remains around 0 only after 10 lags.This suggests that HMC achieves greater independence between samples sooner, indicating potentially more efficient sampling and better exploration of the distribution space compared to the algorithm used in [27].Additionally, an examination of the trace plots in Figure 8 reveals consistent and stable behavior in the MCMC chains, showcasing no discernible oscillation, drift, or other abnormalities.This pattern suggests a thorough exploration of the parameter space, indicating successful convergence to the target distribution.
The small variance observed in Table 1 further supports the reliability of the results.Furthermore, the MCSE is 0 for all βs, indicating that the estimated values are likely close to the true values, as shown in Table 2.The sampler parameters for our model were configured with a n leap f rog value of approximately 50 for all chains, indicating that each iteration of the sampler took around 50 leapfrog steps to explore the parameter space.
From the obtained results, the trace plots do not seem to indicate non-stationary behavior.
Furthermore, upon comparing the results presented in Table 2 for 300 observations with those for 3000 observations, we observe an improvement in estimation accuracy with higher observation counts, indicated by lower RMSE values and reduced variance across most quantiles.This phenomenon is expected, as larger observation sizes generally lead to improved estimations.

Conclusions
This article introduces a Bayesian methodology for simultaneous regression quantile estimation.The framework assumes an ALD for the error term in the quantile regression model and leverages the inter-quantile relationships within the ALD.This approach allows for a comprehensive characterization of the model likelihood across a range of quantiles, providing a nuanced understanding of the estimation process.
Diverging from the methodology presented in [21], which addresses simple quantiles and regression quantiles across diverse data types, our proposed approach is tailored specifically to ALD data in the context of simultaneous regression quantiles.The utilization of HMC for estimation yields promising results, as demonstrated by our RMSE interpretation.However, it is important to note that the current formulation of the HMC method, as described in this study, is not explicitly designed to handle categorical predictors, making it an interesting area for future investigation.
Our estimation procedure, facilitated by HMC, enables the fitting of parametric quantile functions in both univariate linear and multivariate linear quantile settings.The simulation study validates the convergence of HMC algorithms, supported by diagnostic results.Our study demonstrates that HMC achieves greater independence between samples sooner compared to the algorithm used in [27], suggesting potentially more efficient sampling and better exploration of the distribution space.This finding highlights the effectiveness of HMC in Bayesian analysis and underscores its potential for improving computational efficiency in future research.This work represents a substantial advancement towards achieving a fully Bayesian estimation of multiple quantiles.
An exciting perspective to explore for future research involves extending the methodology introduced in [27] to accommodate any conditional distribution of Y|x.Such an extension would broaden the applicability of the proposed approach, enhancing its adaptability to various modeling scenarios.Another perspective within this framework involves conducting sensitivity analyses on key parameters and assumptions related to the HMC algorithm, such as step size and the number of leapfrog steps, to provide further insights into its robustness and effectiveness across various scenarios.It is noteworthy that for both the Metropolis-Hastings and Gibbs algorithms using R, we received the following message: "Cannot compute Gelman & Rubin's diagnostic for any chain segments for variables β 0 .This indicates convergence failure".We received the same message for β 1 .By computing R by hand, we obtained R values less than 1, indicating that the standard deviation between chains is less than the standard deviation within chains.Hence, we did not achieve convergence of the chains.This is not the case when using HMC.Therefore, if we are interested in estimating a univariate linear quantile (not simultaneously), it is preferable to use HMC based on the results obtained in Table A1.Additionally, the estimated quantiles using HMC are better than the ones estimated using Gibbs and Metropolis-Hastings.Hence, for a more complex model like simultaneous regression quantile, it is definitely better to use HMC.This is the main motivation behind using HMC.

Figure 2 .
Figure 2. Histograms showcasing diverse β values for the univariate model.

Figure 4 .
Figure 4. Trace plots depicting various β values for the univariate model.

Figure 6 .
Figure 6.Histograms showcasing diverse β values for the multivariate model with 300 observations.

Table 1 .
Mean, standard deviation, theoretical quantile, and RMSE of β for univariate model with Figure 1.R for univariate model.

Table 2 .
Trace plots depicting various β values for the multivariate model with 300 observations.Mean, standard deviation, theoretical quantile, and RMSE of β for the multivariate model with 300 and 3000 observations.

Table A1 .
Mean, standard variation, the theoretical quantile and RMSE of β for univariate linear model using different algorithms.