Polytomous logistic regression combined with spline smoothing gives a powerful tool for Bayesian density estimation. Using fast array algorithms, multiple dimensions can be handled in a fast and uniform way. The Langevin-Hastings algorithm allows efficient sampling from the associated (re-parameterized) posterior distribution. Illustrations of density estimation are provided, as well as a new approach to smooth quantile regression.
Lambert, P., & Eilers, P. H. C. (2006). Bayesian multi-dimensional density estimation with P-splines (Stat Discussion Paper 0612). https://hdl.handle.net/2078.5/33557