Bayesian hyperparameter learning in a multi-ouput Gaussian Process Regression

Let's imagine I have the following equation $$y_t=f(x_t)+e_t$$, where $$f(x)$$ follows a gaussian process, and $$e_t\sim N(0,\Sigma)$$.

How does one go about to learn the hyperparameters, i.e., $$\Sigma$$ and those of of the gaussian process, in a Bayesian way?

More specifically, how would we write the $$p(Y|X, \theta)$$, where $$\theta$$ are the hyperparameters of the model, and $$Y=(y_1,...,y_N)$$ and $$X=(x_1,...,x_N)$$?

• Could you elaborate a bit more on the MCMC method. How would you write the likelihood? I would like to draw from $p(\theta|X,Y)\propto p(X,Y|\theta)p(\theta)$. My problem is the $p(X,Y|\theta)$, where we've marginalized the function $f$ out... if we had $p(X,f(X), Y|\theta)$ instead, I would know how to find it's expression. Oct 12 '18 at 20:40
• I'm not sure I understand. $Y$ is completely determined by $X$ and $e_t$ in your example. In that example, it's extremely easy to write the likelihood, as X is also an MVN, you'd just be adding $\Sigma$ to the GP's covariance matrix. So, if $e_t \sim N(0, \sigma_n^2 \mathcal I)$, you could simply add the $\sigma_n^2 1(x = x')$ to your kernel Oct 12 '18 at 21:36
• Sorry, I meant $f(X)$ is a multivariate normal, as it's a GP. X is a constant variable right? So, both $f(x_t)$ and $e_t$ are MVNs, hence their sum is as well. Their sum is a GP as well. Oct 12 '18 at 21:48