Suppose I am dead-set on using Bayesian inference on independent and identically distributed data, but I'm lazy and insist on using a parametric likelihood function come what may. I'd be reassured to know that as long as I stick to a particular class of parametric densities I can be sure to recover the true quantiles.
Is there a recipe to construct a parametric family given a vector $q \in (0, 1)^k$ such that the inferred quantiles at $q_1, \dots, q_k$ line up with the true quantiles at those points?
Here's my stab at formalizing this.
To keep things simple I don't mind working with probability distributions over the real line with continuous density with respect to Lebesgue measure, a class we can call $\mathcal{P}$.
The restriction to parametric models can be done by assuming that $\mathcal{G} \subset \mathcal{P}$ can be indexed (smoothly) by a compact subset of $\mathbb{R}^d$, for $d$ finite (so that $\mathcal{G}$ is a smooth manifold as per ArthurB's comment). We can call the parameters $\theta \in \mathbb{R}^d$.
Finally, we have that the Bayesian machinery will converge to the so-called "pseudo-true" parameter values, meaning that we eventually converge to the value $\theta^*$ minimizing the Kullback-Leibler divergence $$\int_{-\infty}^{\infty} \log{\left(\frac{f(x)}{g_{q,\theta}(x)}\right)}f(x) dx$$ for true distribution $F \in \mathcal{P}$ with density $f(x)$.
Now the question reads: for $k > 1$ and for any $F \in \mathcal{P}$ does there exist $\mathcal{G_q} \subset \mathcal{P}$ with $G_{\theta^*}^{-1}(q_j) = F^{-1}(q_j)$ for all $j = 1, \dots, k$ ?
(Here $G^{-1}$ denotes the inverse cumulative distribution function of the distribution $G \in \mathcal{G}$ with density $g(x)$, and similarly for $F$.)
Further, to get to the heart of the issue, I need $\mathcal{G_q}$ to be computable given $q$ as input.
I suspect there is a non-existence result of some kind to be had here. Note that for $k = 1$ this is actually doable, which is the motivation behind the use of the asymmetric Laplace density for this purpose. The trick used there is that the likelihood is maximized at the sample quantile, which is consistent for the true quantile. This does not work naively in the case of $k>1$ because the normalizing constant messes things up so that the parameters are still consistent for the quantiles but they no longer "represent" the corresponding quantiles of the density being used.
