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)).
- cvm is covariance matrix (in the demos always just a scalar>0). It decides how far away the influence of a measurement at x is felt. Small cvm makes the regression "bendier" large cvm "stiffer".
- sigmaf enters mostly as a scaling factor for sigma2 of the measurements. For larger sigmaf the range of linear interpolation is larger. Far outside the x area where measurement were taken, the estimate is mean=0, std=sigmaf.
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');
