Bayesian inference for non-Gaussian errors Following from a previously unanswered question, regression tasks involving measurements with normally distributed noise apply Gaussian processes. But are there any recommended approaches for regression when the uncertainty is not necessarily Gaussian but may follow an arbitrary distribution function? References to relevant literature and especially numerical packages tackling this problem would be greatly appreciated. 
 A: There are several ways in modeling a non-Gaussian error distribution in a Bayesian fashion. But ultimately, it's always how you think a priori about the distribution of your dependent variable/errors.
I would recommend you the following: look at kernel density plots of your dependent variable. Is your dependent variable bounded on a specific interval (e.g. only positive values, etc.pp.)? Is it possible to transform your variable taking logarithm or differences (depending whether you have cross-sectional or time series data)? If you can by any chance transform it to a stochastic process similar to a Normal you will make your life easier.
Otherwise, you may want to look into the BMA (Bayesian Model Averaging) literature. They explicitly take model uncertainty into account, by conditioning on the model they use,
\begin{equation}
p(\theta| y, M_k) = \frac{p(y|\theta, M_k)p(\theta|M_k)}{p(y|M_k)}, \qquad k = 1,...,K
\end{equation}
where the difficulty lies in computing $p(y|M_k)$, the marginal likelihood of the model. Most packages in the realm of BMA are looking at model uncertainty regarding the inclusion/exclusion of exogenous variables, but in principle this can be extended/adapted to specify a set of specific distributions for the errors. Unfortunately, I am not aware of literature specifying different distributional assumptions. A good starting point in the BMA literature may be Hoeting et. al. (1999).
Another approach would be to use mixture models. They specify more than on distribution on the data with weights attached to them. They can be represented in a general way, s. t. $p(y_i) = \sum^K_{k=1} \eta_k p(y_i | \theta_k)$, where
\begin{equation}
y_i = \begin{cases}
T(\theta_1) \text{ if } S_i = 1\\
... \\
T(\theta_K) \text{ if } S_i = K
\end{cases},
\end{equation}
where $T(\cdot)$ is any distribution with a set of parameters $\theta_k$. For example, if you have a multimodal distribution a mixture model with distinctive means may be appropriate. Also distributions with heavy skewness can be characterized by a mixture of normals with distinctive variances. But you can also specify different distributions for each mixture component. A good reference is the monograph by Frühwirth-Schnatter (2006).
I would strongly recommend you to look at your data and start from there in thinking about the statistical assumptions you want to make.
References:
Frühwirth-Schnatter, Sylvia (2006) Finite Mixture and Markov-Switching Models, New York: Springer-Verlag.
Hoeting, J. A., Madigan, D. Raftery, A. E. and Volinsky, C. T. (1999) Bayesian Model Averaging: A Tutorial, Statistical Science, 14(4): 382-401.
