Contents

Quick demo of GP regression/interpolation

GP = Gaussian Process = linear regression (Gptutorial.pdf for formulas and details)

clear
close all
addpath('../utilities/')
rng(0);

Quantity to estimate

Suppose we have a quantity y(x)~N(f(x),sigma2(x)) with mean f(x) and variance sigma2(x). In the tutorial pdf sigma2 is called sigman^2.

x0=0.25;
noiselevel=2e-1;
f=@(x)1./(1+(x-x0).^2);
sigma=@(x)(1+0.2*cos(3*x)+0.2*cos(3*sqrt(2)*x-2))*noiselevel;
xvec=linspace(-1,3,1000);
xrev=xvec(end:-1:1);
clf
lgo=fill([xvec,xrev],[f(xvec)-sigma(xvec),f(xrev)+sigma(xrev)],'c','edgecolor','none');
hold on
grid on
lgo(2)=plot(xvec,f(xvec),'b-');
xlabel('x');
ylabel('f(x)+N(0,yvar(x))');
lgo=lgo([2,1]);
lgt={'mean f','std sigma'};
legend(lgo,lgt,'location','best')
title('Unknown variable y with mean f(x) and std sigma(x)')

Measurements

In measurements we only get an estimate yvar of sigma, say we estimate yvar=noiselevel^2. We make measurements at a few points xp, with estimated value yp and error (std) estimate=noiselevel.

xp=x0+(-0.2:0.4:2);
np=length(xp);
xp=xp+0.2*rand(1,np)-0.1;
yp=f(xp)+sigma(xp).*randn(1,np);
tmp=plot([xp;xp],[yp-noiselevel;yp+noiselevel],'r-');
lgo(end+1)=plot(xp,yp,'ko','linewidth',2);
lgo(end+1)=tmp(1);
lgt=[lgt,{' y est','y std est'}];
legend(lgo,lgt,'location','best')

GP estimate

We assume that the value of y at any x1 is correlated with its value at any x2, and that the correlation is a bell curve of the distance d=|x1-x2| (GPtutorial.pdf formula (1)): k(d)=sigmaf^2*exp(-d^2/(2cvm)).

The assumption uniquely determining the probability distribution for y at a point x, given measurements yp at xp is (3): the joint probability distribution is a multivariate Gaussian with zero mean and covariance matrix given by k(dij) for the points xij. The matrix algebra is implemented in function GPcreate. The interpolation can be evaluated by using GPeval.

yvar=zeros(size(yp))+noiselevel^2;
gp=GPcreate({'cvm',1,'sigmaf',1},xp,yp,yvar);
[y_gp,yv_gp]=GPeval(gp,xvec);
lgo(end+1)=plot(xvec,y_gp,'k','linewidth',2);
tmp=plot(xvec,[y_gp+sqrt(yv_gp);y_gp-sqrt(yv_gp)],'k','linewidth',1);
lgo(end+1)=tmp(1);
lgt=[lgt,{'gp mean','gp std'}];
legend(lgo,lgt,'location','best')

Extrapolation & root/extremum-finding

On the interpolation one can do numerics with full accuracy. For
example, the maximum of the interpolation is

xr=NSolve(@(x)GPevalDiff(gp,x),xp(1),'damping',0.2); % simple Newton iteration, self-written
lgo(end+1)=plot(xr,GPeval(gp,xr),'ks');
lgt=[lgt,{'gp max'}];
legend(lgo,lgt,'location','best')

it=1, |cor|=7e-01, |res|=9e-01, acc=1, renew J=1
it=2, |cor|=4e-01, |res|=7e-01, acc=1, renew J=1
it=3, |cor|=1e-01, |res|=2e-01, acc=1, renew J=1
it=4, |cor|=2e-02, |res|=4e-02, acc=1, renew J=1
it=5, |cor|=2e-03, |res|=4e-03, acc=1, renew J=0
it=6, |cor|=1e-04, |res|=2e-04, acc=1, renew J=0
it=7, |cor|=3e-06, |res|=6e-06, acc=1, renew J=0

Improving the estimate of the maximum

However, the uncertainty is huge. One can use the variance information from the GP regression to look for good x to make additional measurements in. One possible simple criterion is for a new point xn:

xn=xr+linspace(-1,1,30);
[yn,ynvar]=GPeval(gp,xn);
xneffect=NaN(size(xn));
for i=1:length(xn)
    yc=yn(i)-sqrt(ynvar(i));           % new measurement yn+sqrt(ynvar)
    gp_changed=GPcreate(gp,xn(i),yc,noiselevel^2);
    xneffect(i)=GPevalDiff(gp_changed,xr); % changed value of estimated f' at xr by
end                                    % how much?
[~,ix]=sort(abs(xneffect),'descend');
xnew=xn(ix); % in this order the measurements should be taken

Take additional measurements

until the sensitivity of residual in xr to new measurements is smaller than tol

res=@(g,x)GPevalDiff(g,x);
err=FuncSensitivity(gp,xnew,xr,res);
err0=err(1);
k=0;
figure(2);clf
kmax=20;
tol=noiselevel/20;
while k<kmax && err(1)>tol
    k=k+1;
    xreq=xnew(1);
    ynew=f(xreq)+sigma(xreq).*randn(1); % measurement
    gp=GPcreate(gp,xreq,ynew,noiselevel^2);% add new measurement to gp
    xr=NSolve(@(x)GPevalDiff(gp,x),xp(1),'print',0);
    fprintf('new point at x=%g, new estimate for xr=%g\n',xreq,xr);
    GPplot(gp);
    drawnow
    [err,xnew]=FuncSensitivity(gp,xnew,xr,res);
end

new point at x=-0.138109, new estimate for xr=0.487319
new point at x=1.10327, new estimate for xr=0.40078
new point at x=-0.276041, new estimate for xr=0.375693
new point at x=-0.207075, new estimate for xr=0.397944
new point at x=0.965339, new estimate for xr=0.349433
new point at x=0.896373, new estimate for xr=0.35916
new point at x=-0.207075, new estimate for xr=0.325937
new point at x=0.827408, new estimate for xr=0.334859
new point at x=-0.207075, new estimate for xr=0.33105
new point at x=0.827408, new estimate for xr=0.278169
new point at x=-0.345006, new estimate for xr=0.277148
new point at x=0.689477, new estimate for xr=0.299873
new point at x=-0.069144, new estimate for xr=0.299889
new point at x=0.758442, new estimate for xr=0.303351
new point at x=-0.000178439, new estimate for xr=0.298235

Final estimate of maximum

hold on
plot(xr,GPeval(gp,xr),'ks');
xlabel('x');
ylabel('y')
title('approximate maximum of y, interpolated by gaussian Process');

Published with MATLAB® 7.14