How to obtain confidence bounds using nlarx with idGaussianProcess ?

Hi!

I am trying to build a surrogate model with a Non-linear Auto-Regressive with eXogenous inputs model. To do so, I use the system identification toolbox and particularly the matlab functions 'nlarx' using a Gaussian process as output function. I read in the documentation that nlarx with 'idGaussianProcess' uses the matlab function 'fitrgp'. It should be possible than to obtain confidence bounds of the prediction (as it is possible to do so with function '[ypred,ysd,yint] = predict(gprMdl,Xnew)' where gprMdl was obtained with 'fitrgp'). Could you help ?

Here is the code I have actually and I would like to add the confidence bounds of the prediction :

% Load estimation and validation datasets ********************************

z = data_est; % iddata object

z_val = data_val; % iddata object

% Plot estimation and validation datasets *********************************

figure

plot(z)

figure

plot(z_val)

% Options for nlarx command ***********************************************

opt = nlarxOptions;

opt.Display = 'on';

opt.Focus = 'prediction'; % default

% Build nlarx model *******************************************************

% Orders (equivalent to linear regressors)

nbIn = 3;

nbOut = 1;

nanbnk = [1*ones(nbOut,nbOut), 1*ones(nbOut,nbIn), ones(nbOut,nbIn)];

% Initialize idnlarx model and choose the Output function

sysMIMO = idnlarx(z.OutputName, z.InputName, nanbnk, 'idGaussianProcess');

% Estimate nlarx model ****************************************************

sysMIMO = nlarx(z,sysMIMO, opt);

% Compare model response to estimation and validation data ****************

% Estimation data

figure; compare(z,sysMIMO)

% Validation data

figure; compare(z_val,sysMIMO)

[Océane]

I tried something without confidence because I couldn't find any documentation about this use :

[ypred,ysd,yint] = predict(sysMIMO,z_val)

Surprisingly it does ot return any error. But ypred does not correspond to the predicted signal obtained with the 'compare(z_val,sysMIMO)', it does better... Do you have an idea why ? Moreover it seems to work with an output function idLinear, but I can't understand why, it is not a gaussian process. Last thing, I do not understand what is in ysd and yint is returned empty.

Answer 1

[Saarthak Gupta]

Hi Oceane,

It looks like you are trying to obtain the confidence bounds of predictions given by an estimated NARX model with GPR (Gaussian Process Regression) mapping function.

If the underlying ‘RegressionGP’ object of the mapping function was accessible to the user, we could apply the ‘predict’ (compactregressiongp.predict) method of the Statistics and Machine Learning Toolbox to directly compute the standard deviations and 95% prediction intervals of the response variable. Unfortunately, the ‘idGaussianProcess’ mapping function doesn’t allow such access.

However, if we somehow obtain the estimated parameters of the Gaussian Process Regression, we can calculate the probability density of a response y at a point x given the training data as follows:

K (Gaussian kernel), beta, h matrices can be obtained from ‘sys.OutputFcn’ and sigma (noise variance) can be obtained from ‘sys.NoiseVariance’. Here ‘sys’ is the ‘idnlarx’ object returned by the ‘nlarx’ function.

Please note that this expression is valid when the method used for estimating the parameters of the ‘idGaussianProcess’ nonlinear model is ‘exact’ or ‘sd’. Similar expressions can be obtained for ‘sr’ and ‘fic’ fitting methods.

We observe that the responses are normally distributed. Confidence intervals for a particular response can be obtained using standard mathematical techniques for normal distributions. Please refer to the following answer for a detailed explanation: https://www.mathworks.com/matlabcentral/answers/414039

Also, the ‘predict’ method you are using is the one available in the System Identification Toolbox and not the Statistics and Machine Learning Toolbox. They are fundamentally different prediction functions with different inputs and outputs.

The ‘predict’ function in your script returns the estimated values for initial conditions ic and a predictor model sys_pred along with the predicted output yp.

The ‘predict’ function in Statistics and Machine Learning Toolbox: [ypred,ysd,yint] = predict(gprMdl,Xnew) returns the standard deviations ysd and 95% prediction intervals yint of the response variable, evaluated at each observation in Xnew using the trained GPR model

Please refer to the following MATLAB documentation for further reference:

• idGaussianProcess: https://www.mathworks.com/help/ident/ref/idgaussianprocess.html
• predict (System Identification Toolbox): https://www.mathworks.com/help/ident/ref/predict.html
• predict (Statistics and Machine Learning Toolbox): https://www.mathworks.com/help/stats/compactregressiongp.predict.html
• idnlarx: https://www.mathworks.com/help/ident/ref/idnlarx.html
• nlarx: https://www.mathworks.com/help/ident/ref/nlarx.html
• Exact GPR Method: https://www.mathworks.com/help/stats/exact-gpr-method.html
• RegressionGP: https://www.mathworks.com/help/stats/regressiongp.html

Hope this helps!

Best Regards,

Saarthak