Quantile regression is often used when a comprehensive relationship between a response variable and one or more explanatory variables is desired. The traditional frequentists' approach to quantile regression has been well developed around asymptotic theories and efficient algorithms. However, not much work has been published under the Bayesian framework. One challenging problem for Bayesian quantile regression is that the full likelihood has no parametric forms. In this paper, we propose a Bayesian quantile regression method, the linearly interpolated density (LID) method, which uses a linear interpolation of the quantiles to approximate the likelihood. Unlike most of the existing methods that aim at tackling one quantile at a time, our proposed method estimates the joint posterior distribution of multiple quantiles, leading to higher global efficiency for all quantiles of interest. Markov chain Monte Carlo algorithms are developed to carry out the proposed method. We provide convergence results that justify both the algorithmic convergence and statistical approximations to an integrated-likelihood-based posterior. From the simulation results, we verify that LID has a clear advantage over other existing methods in estimating quantities that relate to two or more quantiles.