Estimating correlation hyperparameters of a Gaussian Process

Question

I have an actual function that I need to simulate using a GP model. I've not done this before so I'm unclear of the steps. I have used the true function at different values of the inputs ($\vec X1, \vec X2$) to get some actual data (I think this is my prior distribution).

If I'm not mistaken, the next step is to estimate the parameters. I believe the estimates of the mean and variance are

$$\hat \mu = {{\vec 1\vec R^{-1}\vec y }\over {\vec 1\vec R^{-1}\vec 1}}$$

$$\hat \sigma ^2 = {{(\vec y-\vec 1\hat \mu)^T\vec R^{-1}(\vec y-\vec 1\hat \mu) }\over {n}}$$

Where $y$ is the actual function data, $\vec R$ is the correlation function (which depends on two of the parameters I need to estimate, $\theta_1$ and $\theta_2$).

My understanding is that I need to plug these values into the likelihood or density function and then optimize to get the MLE for $\theta_1$ and $\theta_2$.

I guess this would mean differentiating the likelihood function with respect to $\theta_1$ and $\theta_2$ (yikes) and setting the derivative = 0. I really hope there's an easier way because I don't know how to do that.

I have a feeling I can use the optim() function in R, but I don't know how to use it (I tried unsuccessfully).

Can anyone help?

I think once I have these estimates I can stick them into the conditional mean and variance formulas and go from there.

Edit: An alternative I think is to use MCMC sampling to estimate the integral of the posterior. Is that the way to go? If so are there any good tutorials on how to actually do it?

Answer 1

Ok first off what you are doing is akin to the backbone of bayesian optimisation and quadrature. That is trying to learn a true noise free function from observations.

However you seem quite confused and I think it's worth pointing you to a good reference rather than writing a massive post.

Implement the GP: Look at sudo code page 19 of this

Optimise Hyperparameters: Equation 5.8 and 5.9 of this chapter shows you the log marginal likelihood and its derivatives with respect to parameters. You need to maximise this (or minimise the negative of it). Depending on the kernel you wish to use the derivatives vary in their difficulty. However I show how to optimise the RBF kernel in this answer.

To be honest I would not use MCMC to optimise hyper parameters, especially if this is your first GP. Instead use a convex optimisation algorithm as there a many off the shelf blackboxes that work fine for simple GPs.

I can't be much help re: R. My advice though is to look at GPy, GPML or even pyGPs. Read through there code for the basic samples and read through the code of the library to see the flow of execution and how everything fits together. It's a massive help!