A black box function $f: \mathbb{R}^n \rightarrow \mathbb{R}$, which is evaluated pointwise subject to Gaussian noise, i.e., $f(x) + \mathcal{N}(\mu(x),\sigma(x)^2)$, can be minimized using Bayesian optimization where a Gaussian Process is used as a noisy function model.
How can Bayesian optimization be used for functions subject to non-Gaussian noise, e.g., skewed distributions?
Are there any implementations that support this setting?
There are Gaussian process models with non-Gaussian likelihood: The prior distribution on the function $f$ is still a Gaussian process but the noise term is not Gaussian anymore, i.e. the likelihood $p(y | f)$ is not assumed to be Gaussian anymore. As a consequence the analytical results are lost and drawing inference now requires approximation methods such as MCMC or Laplace approximation.
For several distributions this is implemented and explained as part of the GPML Matlab package, available and explained here. The table of inference methods in section 3d ("A More Detailed Overview") gives an overview of what distributions have been implemented for the likelihood and what inference method is available for each of them.
The only articles that I can link you to right now (because I bookmarked them at some point) are on the Student's $t$ distribution:
