  In Bayesian applications, there is a huge interest in rapid and accurate
estimation of the posterior distribution, particularly for high dimensional or
hierarchical models. In this article, we propose to use optimization to solve
for a joint distribution (random transport plan) between two random variables,
$\theta$ from the posterior distribution and $\beta$ from the simple
multivariate uniform. Specifically, we obtain an approximate estimate of the
conditional distribution $\Pi(\beta\mid \theta)$ as an infinite mixture of
simple location-scale changes; applying the Bayes' theorem,
$\Pi(\theta\mid\beta)$ can be sampled as one of the reversed transforms from
the uniform, with the weight proportional to the posterior density/mass
function. This produces independent random samples with high approximation
accuracy, as well as nice theoretic guarantees. Our method shows compelling
advantages in performance and accuracy, compared to the state-of-the-art Markov
chain Monte Carlo and approximations such as variational Bayes and normalizing
flow. We illustrate this approach via several challenging applications, such as
sampling from multi-modal distribution, estimating sparse signals in high
dimension, and soft-thresholding of a graph with a prior on the degrees.
