How to use the gaussian process regression function in matlab 2015b ?

I'm trying to reproduce the picture at the end of the documentation: http://www.mathworks.com/help/stats/gaussian-process-regression-models.html the one on the left.

However, I cant quite do it.

This is what I have so far:

fun = @(x) x.*sin(x);

xx = linspace(0,10,100)';

yy = fun(xx);

nd = 5;
xd = linspace(1,9,nd)';
yd = fun(xd);

gp = fitrgp(xd,yd,'KernelFunction','squaredexponential');

figure(1); clf
plot(xx,yy, 'g-')
hold on
plot(xx,predict(gp,xx), 'r-')

However this doesn't work. It looks like the beta parameter is set at a constant and you can change the basis function used. But, why use Gaussian Processes if you have to provide it with the function you're trying to emulate? I'm trying to use GPs to model simulation data and the process that generate them can't be written as a nice function (basis function).

Answer 1

[Gautam Pendse]

Hi Regis,

Here's how you can reproduce the doc example:

      % True curve.
      fun = @(x) x.*sin(x);
      xx  = linspace(0,10,100)';
      yy  = fun(xx);

      % Example points (from doc). You can do xd = linspace(1,9,nd)' as well 
      % but you will get a different fit.    
      xd = [1,3,5,6,7,8]';
      yd = fun(xd);

      % Fit a GP model. Initialize 'Sigma' to a small value. A GP estimates
      % its parameters by maximizing the marginal log likelihood. Depending
      % on the data, the marginal log likelihood can have multiple local
      % optima corresponding to different interpretations of the data.
      % Initializing 'Sigma' to a small value discourages the high noise
      % variance interpretation of the data.
      gp = fitrgp(xd,yd,'KernelFunction','squaredexponential','sigma',0.1,'verbose',1);

      % Plot.
      figure(1); clf
      plot(xx,yy, 'r-.')
      hold on;

      [ypred,~,yint] = predict(gp,xx);

      plot(xx,ypred, 'g-');
      plot(xx,yint(:,1),'k-');
      plot(xx,yint(:,2),'m-');

      plot(xd,yd,'ro');

      legend('f(x) = x.*sin(x)','GPR predictions','Lower 95% interval','Upper 95% interval','Observations','Location','Best');

Hope this helps,

Gautam

[Regis]

Thank you very much for your answer. It works. Quick follow up, and I'm sorry if this is obvious but I cant figure this out.

How do you specify the noise in the data? i.e and make the figure on the right-hand side of that documentation? I figured out that doing something like yd = fun(xd) + normrnd(0,1,6,1); works and give me a picture similar to the one on the left but not always (depending on the normrnd I get). How do I set the noise for each data point? The theory says it's a vector that you add to the diagonal of the covariance matrix of the GP, but when I try to set that with the kernelParameter variable, it doesn't work. Thanks again for your help

Answer 2

[Gautam Pendse]

Replace xd and yd by something like this:

xd = linspace(0,10,20)';
yd = fun(xd) + 1*randn(20,1);

Rest of the code remains the same.

Gautam