Download the housing data [1], from the UCI Machine Learning Repository [4].

The dataset has 506 observations. The first 13 columns contain the predictor values and the last column contains the response values. The goal is to predict the median value of owner-occupied homes in suburban Boston as a function of 13 predictors.

Load the data and define the response vector and the predictor matrix.

load('housing.data');
X = housing(:,1:13);
y = housing(:,end);

Fit a GPR model using the squared exponential kernel function with separate length scale for each predictor. Standardize the predictor variables.

gprMdl = fitrgp(X,y,'KernelFunction','ardsquaredexponential','Standardize',1);

Create a cross-validation partition for data using predictor 4 as a grouping variable.

rng('default') % For reproducibility
cvp = cvpartition(X(:,4),'kfold',10);

Create a 10-fold cross-validated model using the partitioned data in cvp.

cvgprMdl = crossval(gprMdl,'CVPartition',cvp);

Compute the regression loss for in-fold observations using models trained on out-of-fold observations.

L = kfoldLoss(cvgprMdl)

L =

    9.5299

Predict the response for in-fold observations, i.e. observations not used for training.

ypred = kfoldPredict(cvgprMdl);

For every fold, kfoldPredict predicts responses for observations in that fold using the models trained on out-of-fold observations.

Plot the actual responses and prediction data.

plot(y,'r.');
hold on;
plot(ypred,'b--.');
axis([0 510 -15 65]);
legend('True response','GPR prediction','Location','Best');
hold off;

Download the abalone data [2], [3], from the UCI Machine Learning Repository [4] and save it in your current directory with the name abalone.data.

Read the data into a table.

tbl = readtable('abalone.data','Filetype','text','ReadVariableNames',false);

The dataset has 4177 observations. The goal is to predict the age of abalone from 8 physical measurements.

Fit a GPR model using the subset of regressors (sr) method for parameter estimation and fully independent conditional (fic) method for prediction. Standardize the predictors and use a squared exponential kernel function with a separate length scale for each predictor.

gprMdl = fitrgp(tbl,tbl(:,end),'KernelFunction','ardsquaredexponential',...
      'FitMethod','sr','PredictMethod','fic','Standardize',1);

Cross-validate the model using 4-fold cross validation. This partitions the data into 4 sets. For each set, fitrgp uses that set (25% of the data) as the test data, and trains the model on the remaining 3 sets (75% of the data).

rng('default') % For reproducibility
cvgprMdl = crossval(gprMdl,'KFold',4);

Compute the loss over individual folds.

L = kfoldLoss(cvgprMdl,'mode','individual')

L =

    4.3669
    4.6896
    4.0565
    4.3162

Compute the average cross-validated loss on over all folds. The default is the mean squared error.

L2 = kfoldLoss(cvgprMdl)

L2 =

    4.3573

This is equal to the mean loss over individual folds.

mse = mean(L)

mse =

    4.3573