The "Sigma" parameter in the "fitrgp" function is designed for homoscedastic noise, meaning it cannot directly take a vector of output standard deviations (std_output) that varies across data points. To incorporate both the mean and standard deviation when training a surrogate model with "fitrgp", you can use two separate Gaussian Process (GP) models: one to model the mean and another to model the variance. After fitting these models, you can define a custom GP that integrates the variance model, and then fit a final GP using this custom kernel. This methodology is based on the research paper "Improved Most Likely Heteroscedastic Gaussian Process Regression via Bayesian Residual Moment Estimator," which can be found in the following link: https://ieeexplore.ieee.org/document/9103623
Refer to the following steps for a detailed understanding of the approach:
- Fit a GP surrogate model to your mean values as usual:
gprMean = fitrgp(x, mean_output, 'KernelFunction', 'squaredexponential');
- Fit a separate GP to the sample variances (i.e., squared standard deviations)
variance_output = std_output.^2;
gprVar = fitrgp(x, variance_output, 'KernelFunction', 'squaredexponential');
- Create a kernel function that, for the training points, adds the predicted variance from gprVar to the diagonal of the kernel matrix:
kfcn = @(XN, XM, theta) myHeteroKernel(XN, XM, theta, gprVar);
function K = myHeteroKernel(XN, XM, theta, gprVar)
sigmaL = exp(theta(1));
sigmaF = exp(theta(2));
D2 = pdist2(XN, XM, 'squaredeuclidean');
K = (sigmaF^2) * exp(-0.5 * D2 / sigmaL^2);
if isequal(XN, XM) % Only add variance for training data
varXN = predict(gprVar, XN);
varXN = max(varXN, 0);
K = K + diag(varXN);
end
end
- Now, fit the GP to your mean data using the custom kernel:
gprHetero = fitrgp(x, mean_output, ...
'KernelFunction', kfcn, ...
'KernelParameters', [0; 0], ...
'Sigma', 0); % Set Sigma=0, as noise is handled by the kernel
- Make predictions for the input values using:
xq = (5:0.1:13)';
[yhat, ysd, yci] = predict(gprHetero, xq);
As the number of input data points is only 5, the variance GP may be rough. Implementing the above script results the below image:
Refer to the following documentations to know more about "fitrgp" function: https://www.mathworks.com/help/stats/fitrgp.html
