  A broad class of models that routinely appear in several fields can be
expressed as partially or fully discretized Gaussian linear regressions.
Besides including basic Gaussian response settings, this class also encompasses
probit, multinomial probit and tobit regression, among others, thereby yielding
to one of the most widely-implemented families of models in applications. The
relevance of such representations has stimulated decades of research in the
Bayesian field, mostly motivated by the fact that, unlike for Gaussian linear
regression, the posterior distribution induced by such models does not seem to
belong to a known class, under the commonly-assumed Gaussian priors for the
coefficients. This has motivated several solutions for posterior inference
relying on sampling-based strategies or on deterministic approximations that,
however, still experience computational and accuracy issues, especially in high
dimensions. The scope of this article is to review, unify and extend recent
advances in Bayesian inference and computation for this class of models. To
address such a goal, we prove that the likelihoods induced by these
formulations share a common analytical structure that implies conjugacy with a
broad class of distributions, namely the unified skew-normals (SUN), that
generalize Gaussians to skewed contexts. This result unifies and extends recent
conjugacy properties for specific models within the class analyzed, and opens
avenues for improved posterior inference, under a broader class of formulations
and priors, via novel closed-form expressions, i.i.d. samplers from the exact
SUN posteriors, and more accurate and scalable approximations from VB and EP.
Such advantages are illustrated in simulations and are expected to facilitate
the routine-use of these core Bayesian models, while providing a novel
framework to study theoretical properties and develop future extensions.
