Predictions using Hyperparameter Posterior Samples

[mh510]

I encountered an error when I ran the following example code. The code basically computes two GPs, one uses maximum likelihood estimation for the underlying hyperparameters and the other one uses samples of the hyperparameter posterior. It also visualizes the predictive distributions for both GPs.

%% code
clear; close all; clc; rng(42);

% - define ground truth
xlb = 0.5;
xub = 2.5;
truefcn = @(x) sin(10*pi*x)./(2*x)+(x-1).^4;      % gramacy and lee testfunction
xgrid = linspace(xlb, xub, 1000)';
fgrid = truefcn(xgrid);

% - create data
N = 444;
s2n = 0.01;
x = sort(xlb+(xub-xlb)*rand(N, 1));
y = truefcn(x)+sqrt(s2n)*randn(N, 1);

% - use maximum likelihood
gp1 = gp_set('lik', lik_gaussian('sigma2', 0.1), 'cf', gpcf_sexp('lengthScale', 1, 'magnSigma2', 0.1));
gp1 = gp_optim(gp1, x, y);
[~, ~, ~, ymean1, yvar1] = gp_pred(gp1, x, y, xgrid);

% - use posterior sampling
lik = lik_gaussian('sigma2', 1, 'sigma2_prior', prior_gaussian('mu', 0.01, 's2', 1));
gpcf = gpcf_sexp('lengthScale', 1, 'magnSigma2', 1, ...
                 'lengthScale_prior', prior_gaussian('mu', 1, 's2', 1), 'magnSigma2_prior', prior_gaussian('mu', 0.01, 's2', 1));
gp2 = gp_set('lik', lik, 'cf', gpcf);
gp2 = thin(gp_mc(gp2, x, y, 'nsamples', 25), 5);
[~, ~, ~, ymean2, yvar2] = gp_pred(gp2, x, y, xgrid);            % <--- triggers the error

% - visualize
figure(1);
hold on; grid on;
plot(xgrid, fgrid, 'g');
plot(x, y, 'k+');
plot(xgrid, ymean1, 'b');
plot(xgrid, ymean2, 'r');
plot(xgrid, ymean1+3*sqrt(yvar1), 'b--');
plot(xgrid, ymean1-3*sqrt(yvar1), 'b--');
plot(xgrid, ymean2+3*sqrt(yvar2), 'r--');
plot(xgrid, ymean2-3*sqrt(yvar2), 'r--');

After debugging I found the corresponding error-message in gpmc_preds.m:

if nargout > 2 && isempty(yt)
   error('mc_pred -> If lpy is wanted you must provide the vector yt as an optional input.')
end

My workaround to make things work is to uncomment this if-condition and adjust the code around line 238 like this:

if ~isempty(yt)
   [Ef(:,i1), Varf(:,:,i1), lpy(:,i1), Ey(:,i1), Vary(:,:,i1)] = gp_pred(Gp, x, y, xt, options);
   else
   [Ef(:,i1), Varf(:,:,i1), ~, Ey(:,i1), Vary(:,:,i1)] = gp_pred(Gp, x, y, xt, options);  
   lpy=[];
end

[jpvanhat]

Thank you for pointing this out. I just made a pull request to fix this for "develop" branch.