Matlab GPR -> equation for prediction

Dear MathWorks-Team

I have a question regarding a GPR model (classreg.learning.regr.CompactRegressionGP). I trained a model in Matlab which works very well. The problem is, I want to use it in javascript for a webpage and this is where my issue starts. I know, that there is possibility to export a model using Matlab Coder into c++/c this does not help in my case since I really need to implement it myself in the context of this webpage.

My knowledge on GP is however only basic and I struggle to understand the theory here (https://ch.mathworks.com/help/stats/exact-gpr-method.html). I know that I can derive the hyperparameters, ect from here (https://ch.mathworks.com/help/stats/classreg.learning.regr.compactregressiongp-class.html). The issue is, this is for someone who is not very very deep in the material impossible to understand/implent himself.

Could you therefore please write me down an easy equation based on the parameters of a model, how I can rewrite the prediction function?

This is my model:

obj_TrainedModel = 

  classreg.learning.regr.CompactRegressionGP
           PredictorNames: {'column_2'  'column_3'  'column_4'}
             ResponseName: 'Y'
    CategoricalPredictors: []
        ResponseTransform: 'none'
           KernelFunction: 'Matern52'
        KernelInformation: [1×1 struct]
            BasisFunction: 'Constant'
                     Beta: 0.5967
                    Sigma: 0.2731
        PredictorLocation: [3×1 double]
           PredictorScale: [3×1 double]
                    Alpha: [39×1 double]
         ActiveSetVectors: [39×3 double]
            PredictMethod: 'Exact'
            ActiveSetSize: 39

I appreciate all help a lot! Thanks in advance,

Best regards, J

[Jonas Widmer]

Noone?

Answer 1

[Jonas Widmer]

Found a solution myself. Here is the code if someone else needs it:

clc
load('my_GPR')


RegObj = obj_TrainedSet1.RegressionGP;
SigmaL = RegObj.KernelInformation.KernelParameters(1);  %length scale
SigmaF = RegObj.KernelInformation.KernelParameters(2);  %kernel scale
Sigma = RegObj.Sigma; %measurement noise

Alpha = RegObj.Alpha; % precomputed (K + Sigma*I)^-1*Y

Xraw = table2array(RegObj.X);
Yraw = RegObj.Y;
offsetX = RegObj.PredictorLocation'; % X offset
scaleX = RegObj.PredictorScale';     % X scaling
offsetY = RegObj.Beta;                  % Y offset

% after scaling (this is what goes into the GP)
X = repmat(scaleX,size(Xraw,1),1).\(Xraw - offsetX);
Y = Yraw - offsetY;

dist = @ (x1,x2) sqrt(sum((x1-x2).*(x1-x2),2)); %euclidean distance sqrt(x1-x2)'*(x1-x2)
matern52 = @(x1,x2) SigmaF^2*(1+ sqrt(5)*dist(x1,x2)/SigmaL + 5*dist(x1,x2).^2/(3*SigmaL^2)).*exp(-sqrt(5)*dist(x1,x2)/SigmaL);

K = [];
for i=1:size(X,1)
    K(:,i) = matern52(X(i,:),X); % K matrix
end

Alpha_reconstructed = (K + Sigma^2*eye(size(X,1)))\Y; % this is already computed in Alpha
norm(Alpha-Alpha_reconstructed);

% test prediction of training data
for i=1:size(X,1)
   pred_test(i) = matern52(X(i,:),X)'*Alpha;
end

% hold off
% plot(pred_test)
% hold on
% plot(Y) % fits the training data quite well

% for a general test point
x_new_raw = rand(1,3);
x_new = scaleX.*(x_new_raw - offsetX);
x_new = scaleX.\(rand_Input - offsetX);

mean_prediction = matern52(x_new,X)'*Alpha + offsetY;
var_prediction = matern52(x_new,x_new) + matern52(x_new,X)'*((K + Sigma*eye(size(X,1)))\matern52(x_new,X));