MATLAB: Use of fitrgp to solve simple example (Rasmussen book)

fitrgpgaussian processesgpmlStatistics and Machine Learning Toolbox

I am learning Gaussian Processes and I am trying to replicate the canonical example of Rasmussen book. I understand it should be a simple application of fitrgp, but I cannot get it. I always get a flat response:

My idea is to have a very simple training set in a 2-D example

X=[1,3,7,10]';% those are the given points

y=[17,14,10,21]';% those are the expected responses

%Then I train a model

tbl=table(X,y);

grp=fitrgp(tbl,'y');

%and try to evaluate the fit for the interval [0,20]

xnew=linspace(0,20)';

[ypred,ysd]=predict(grp,xnew);

plot(xnew,ypred, xnew, ypred+ysd)

I was looking for something similar to Rasmussen result

But I only get a flat response

I have been looking in the web and I have found a example code that does the trick, but I understood that fitrgp should include and supersede it.

Am I correct?. Should I use some parameters with fitrgp different from the default? (I was reading/trying some of them carefully and I could not find the right combination)

Thank you in advance

Best Answer

Hi Jose,

Function fitrgp fits the model by maximizing the marginal log likelihood. In some cases, there could be multiple local optima corresponding to different interpretations of the data. Consider initializing the noise standard deviation to a smaller value using the 'Sigma' name/value pair like this:

x = [-4;-2;-1;0;2];
y = [-2;0;1;2;-1];
gpr = fitrgp(x,y,'Sigma',0.1);
plot(x,y,'b+','DisplayName','Data');
hold on;
xtest = linspace(-5,5,1000)';
[pred,~,ci] = predict(gpr,xtest);
plot(xtest,pred,'r','DisplayName','Prediction');
hold on;
plot(xtest,ci(:,1),'c','DisplayName','Lower 95% Limit');
plot(xtest,ci(:,2),'k','DisplayName','Upper 95% Limit');
legend('show','Location','Best');
shg;

Also, have a look at this doc example which illustrates the effect of specifying initial parameter values to fitrgp:

https://www.mathworks.com/help/stats/fitrgp.html#buv_uxs-4

Hope this helps,

Gautam

Related Solutions

MATLAB: How to choose the best kernel for a Gaussian process regression, possibly using “bayesopt” function? How to predict unseen data with the resulting model

There are two ways you could do this. One of them involves "bayesopt" function.

(1) The easiest way is to use the "fitrgp" function. There is an example of how to optimize a fitrgp hyperparameter (Sigma) on the "fitrgp" function Documentation page:

https://www.mathworks.com/help/stats/fitrgp.html?searchHighlight=fitrgp

To optimize ‘KernelFunction’, using default values for all other hyperparameters, you could do the following (using data from the documentation example):

load(fullfile(matlabroot,'examples','stats','gprdata2.mat'))
model1 = fitrgp(x,y,'OptimizeHyperparameters',{'KernelFunction'});

Because fitrgp tends to be sensitive to the Sigma parameter, you will probably get a better result if you simultaneously optimize KernelFunction and Sigma. Also, searching categorical variables inherently requires more function evaluations, so you should probably run it longer:

model2 = fitrgp(x,y,'OptimizeHyperparameters',{'KernelFunction','Sigma'},...
    'HyperparameterOptimizationOptions', struct('MaxObjectiveEvaluations',100));

The "fitrgp" function documentation page lists the hyperparameters that are eligible for optimization, and what values are searched:

https://www.mathworks.com/help/stats/fitrgp.html#input_argument_namevalue_OptimizeHyperparameters.

You can use this model to make predictions on unseen data directly with the "predict" function:

model2 = fitrgp(x,y,'OptimizeHyperparameters',{'KernelFunction','Sigma'},...
    'HyperparameterOptimizationOptions', struct('MaxObjectiveEvaluations',100));
% Predict on unseen data

xNew = 0.2
ypred = predict(model2, xNew)

(2) A second approach would be to use "bayesopt" function directly, defining your own objective function. Here, the objective function is the 5-fold crossvalidation loss of the model.

load(fullfile(matlabroot,'examples','stats','gprdata2.mat')) 
% Define kernel variable
kernel = optimizableVariable('KernelFunction',{'exponential','squaredexponential','matern32','matern52',...
'rationalquadratic','ardexponential','ardsquaredexponential','ardmatern32','ardmatern52','ardrationalquadratic'},...
'Type','categorical') 
% Call bayesopt, capturing x and y in the objective function

bo = bayesopt(@(T)objFcn(T,x,y), kernel)
% Define an objective function
function Loss = objFcn(Vars, x, y)
m = fitrgp(x, y, 'KernelFunction', char(Vars.KernelFunction), 'KFold', 5);Loss = kfoldLoss(m);
end

(2.1) To optimize both KernelFunction and Sigma using "bayesopt" function, you can do the following:

% Define 2 variables
kernel = optimizableVariable('KernelFunction',{'exponential','squaredexponential','matern32','matern52',...
'rationalquadratic','ardexponential','ardsquaredexponential','ardmatern32','ardmatern52','ardrationalquadratic'},...
'Type','categorical')
sigma = optimizableVariable('Sigma',[1e-4,10],'Transform','log')
 
% Call bayesopt, capturing x and y in the objective function
bo = bayesopt(@(T)objFcn(T,x,y), [sigma, kernel], 'MaxObjectiveEvaluations', 100)
 
function Loss = objFcn(Vars, x, y)
m = fitrgp(x, y, 'KernelFunction', char(Vars.KernelFunction), ...
                 'Sigma', Vars.Sigma, 'ConstantSigma', true,...
                 'KFold', 5);
Loss = kfoldLoss(m);
end

To get a model that can make predictions, you need to fit a model without passing the ‘kfold’ argument, but keeping the optimal kernel function obtained from "bayesopt" function.

kernel= bo.XAtMinObjective.KernelFunction
model = fitrgp(x, y, 'KernelFunction', char(kernel));
% Predict on unseen data
xNew = 0.2
ypred = predict(model2, xNew)

MATLAB: Sampling from Posterior Distribution of GPR Model from fitrgp()

Hi Sterling,

Here's an example illustrating how to sample from the posterior distribution of a GPR model. The code uses an undocumented function predictExactWithCov. If you have categorical predictors, you would need to convert them to dummy variables before using that function.