# Estimating correlation hyperparameters of a Gaussian Process

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?