% Example 7.1
x=10*rand(1,10);
xnear=x+0.1; xfar=x+1;
y=sin(x);
ynear=sin(xnear)
yfar=sin(xfar)
rnear=corrcoef(y, ynear) %rnear=0.9894;
rfar=corrcoef(y, yfar) %rfar=0.4229
%-----------------------------------------------------------------------------
%
% Example 7.2
x=linspace(0,4,5);
[col,raw] = meshgrid(x,x);
dist=abs(col-raw);
theta1=2.0;
R1=exp(-(dist/theta1).^2)
det(R1)
theta2=0.1;
R2=exp(-(dist/theta2).^2)
det(R2)
%-----------------------------------------------------------------------------
%
% Example 7.5
y=[0 1 0]';
X=[1 0;1 1;1 2];
R=[ 1 exp(-1) exp(-4)
    exp(-1) 1 exp(-1)
    exp(-4) exp(-1) 1];
Rinv=inv(R);
A=X'*Rinv*X;
Ainv=inv(A);
for i=1:51
 x(i)=0.04*(i-1);
 r = [exp(-((0-x(i))/th)^2);exp(-((1-x(i))/th)^2);exp(-((2-x(i))/th)^2)];
 xi= [1; x(i)];
 w(i,:)=r'*Rinv + (xi'-r'*Rinv*X)*Ainv*X'*Rinv;
end
plot(x,w(:,1),x,w(:,2),x,w(:,3))
legend('w_1','w_2','w_3')
%-----------------------------------------------------------------------------
%
% Example 7.7
y=[0 1 0]';
X=[1 1 1]';
for k=1:100
th(k)=0.02*k;
e1=exp(-(1/th(k))^2); e2=exp(-(2/th(k))^2);
R=[1 e1 e2; e1 1 e1; e2 e1 1];
Rinv=inv(R);
detR=det(R);
beta=(X'*Rinv*y)/(X'*Rinv*X);
ym=y-X*beta;
Like(k)=log((ym'*Rinv*ym/2)^4*detR);
end
plot(th,Like)
%-----------------------------------------------------------------------------
%
% Example 7.8
x=[0 0.5 1 1.5 2];
y=[0 0.7071 1 0.7071 0]';
X=[1 1 1 1 1]';
ny=size(x,2);
R=zeros(ny,ny);
for k=1:100
 th(k)=0.02*k;
 for i=1:ny
  for j=1:ny
   R(i,j)=exp(-((x(i)-x(j))/th(k))^2);
  end
 end
 Rinv=inv(R);
 detR=det(R);
 beta(k)=(X'*Rinv*y)/(X'*Rinv*X);
 ym=y-X*beta(k);
 Like(k)=log((ym'*Rinv*ym/4)^8*detR);
end
figure; plot(th,Like)
%-----------------------------------------------------------------------------
%
% Example 7.9
x=[0 5 10 15 20]';                                            % input variable
y=[1 0.99 0.99 0.94 0.95]';                                          % samples
X=ones(5,1);                                                   % design matrix
ny=length(y); np=size(X,2);
th=5.2;
for k=1:ny; for l=1:ny;
 R(k,l)=exp(-(norm(x(k,:)-x(l,:))/th)^2);
end; end;
Rinv=inv(R);
betaH=(X'*Rinv*X)\(X'* Rinv*y);
sigmaH=sqrt(1/(ny-np)*((y-X*betaH)'*Rinv*(y-X*betaH)));
th=zeros(20,1); Obj=zeros(20,1);
for i=1:20
 th(i)=0.5*i;
 for k=1:ny; for l=1:ny;
  R(k,l)=exp(-(norm(x(k,:)-x(l,:))/th(i))^2);
 end; end;
 Rinv=inv(R);
 betaH=(X'*Rinv*X)\(X'* Rinv*y);
 sigmaH=sqrt(1/(ny-np)*((y-X*betaH)'*Rinv*(y-X*betaH)));
 Obj(i)=log(sigmaH^(2*(ny-np))*det(R));
end
plot(th,Obj,'linewidth',2); grid on;
xNew=10; %or xNew=14
for k=1:ny; r(k,1)=exp(-(norm(x(k,:)-xNew)/th)^2); end;
gpDepar=r'*Rinv*(y-X*betaH);
%-----------------------------------------------------------------------------
%
% Example 7.10
clear; ny=20; nth=10;
x=linspace(0,2*pi,ny)'; thmin=0.1;
y=sin(x)+0.2*x;
th=linspace(thmin,2.0,nth)';
s=[-pi:0.05:3*pi]';
X=ones(ny,1);
for k=1:nth                        % Loop for increasing hyperparameter
 for i=1:ny, for j=1:ny                            % Correlation matrix
  R(i,j)=exp(-((x(i)-x(j))/th(k)).^2);
 end, end
 if rcond(R) < eps, break; end         % Check for reciprocal cond. No.
 Rinv=inv(R);
 beta=inv(X'*Rinv*X)*(X'*Rinv*y);             % Global coefficient beta
 ym=y-beta;
 Like(k)=log((ym'*Rinv*ym/(ny-1))^(ny-1)*det(R)); % Neg. log-likelihood
 for i=1:ny
  r(:,i)=exp(-((s-x(i))/th(k)).^2);   % Correlation at prediction point
 end
 yhat(:,k)=beta+r*Rinv*ym;                        % Kriging predictions
 figure(1); plot(s,yhat); hold on;
 fprintf("k=%2d th=%5.3f beta=%5.3f Like=%8.3e\n",k,th(k),beta,Like(k))
end
plot(x,y,'*',s,sin(s)+0.2*s,'k--'); hold off;
figure(2); plot(th(1:size(Like,2)),Like','b*-');
%-----------------------------------------------------------------------------
%
% Example 7.11
xi=1;
w=Rinv*r+Rinv*X*((X'*Rinv*X)\(xi'-X'*Rinv*r));
ySigmaH=sigmaH*sqrt(w'*R*w-2*w'*r+1);
% using the inverse calculation
gpMean=0.9482; % from Example 7.9
PI=[ gpMean + tinv(0.05,ny-np)*ySigmaH, ...
     gpMean + tinv(0.95,ny-np)*ySigmaH]
% using the random samples
ns=5e3;
tDist=trnd(ny-np,1,ns);
yHat=gpMean+tDist*ySigmaH;
PI=prctile(yHat,[5 95])
%-----------------------------------------------------------------------------
%
% Example 7.12
ny=10; nth=4;
x=linspace(0,2*pi,ny)'; s=[-pi:0.05:3*pi]';
y=sin(x)+0.2*x;
th=[0.1; 0.5; 4.36; 5.5];
X=ones(ny,1);
for k=1:nth
 for i=1:ny, for j=1:ny
  R(i,j)=exp(-((x(i)-x(j))/th(k)).^2);
 end, end
 if rcond(R) < eps, break; end
 Rinv=inv(R);
 beta=inv(X'*Rinv*X)*(X'*Rinv*y);
 sigmaH=ym'*Rinv*ym/(ny-1);
 ym=y-beta;
 for i=1:ny
  r(:,i)=exp(-((s-x(i))/th(k)).^2);
 end
 yhat(:,k)=beta+r*Rinv*ym;
 A=X'*Rinv*X; Ainv=inv(A);
 for i=1:size(s,1)
  w=r(i,:)*Rinv + (1-r(i,:)*Rinv*X)*Ainv*X'*Rinv;
  ySigmaH(i,k)=sigmaH*sqrt(w*R*w'-2*w*r(i,:)'+1);
 end
 figure; hold on; 
 plot(s,yhat(:,k),x,y,'*',s,sin(s)+0.2*s,'r-','LineWidth',1);
 patch([s;flipud(s)],[yhat(:,k)+tinv(0.05,ny-1)*ySigmaH(:,k); ...
       flipud(yhat(:,k)+tinv(0.95,ny-1)*ySigmaH(:,k))], ... 
       'k','FaceAlpha',0.1);
 fprintf("k=%2d th=%5.3f beta=%5.3f Like=%8.3e, detR=%8.3e\n",... 
          k,th(k),beta,Like(k),det(R))
end
%-----------------------------------------------------------------------------
%
% Example 7.13
myfun=@(x) 0.5*((6*x-2).^2.*sin(12*x-4))+10*(x-0.5)+5;   % In-line function
xs=[0.1548,0.7148,0.4934,0.9932,0.5862,0.0007,0.3638,0.8198]';    % Samples
ys=myfun(xs);                                  % Samples from true function
trPRS=0;                                                 % Ordinary Kriging
xp=(0:0.01:1)';                                         % Prediction points
[yp, ysig, nb]=Kriging(xs,ys,xp,trPRS,0);         % Build Kriging surrogate
ny=size(ys,1);                                          % Number of samples
cilb = yp + tinv(0.025,ny-nb)*ysig;  % Lower-bounds of confidence intervals
ciub = yp + tinv(0.975,ny-nb)*ysig;  % Upper-bounds of confidence intervals
ytrue=myfun(xp);                % True function values at prediction points
figure; hold on;                                         % Plotting results
xShaded=[xp;sort(xp,'descend')]; yShaded=[cilb;ciub(end:-1:1)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(xp,yp,'r-',xs,ys,'k+',xp,ytrue,'k-','Markersize',10,'Linewidth',1.5)
xlabel('x');ylabel('y(x)');
legend('95% PI','Kriging prediction','Samples','True function');
%-----------------------------------------------------------------------------
%
% Example 7.14
BraninHu = @(x) (x(:,2)-1.275*(x(:,1)/pi).^2+5*(x(:,1)/pi)-6).^2 ...
                 +10*(1-1/(8*pi))*cos(x(:,1))+10;
rng('default'); N=31; dN=15/(N-1); ny=20;
xs=lhsdesign(ny,2,'iteration',5000);                       % 20 LHS samples
xs=xs.*[15, 15]+[-5, 0];
ys=BraninHu(xs);                               % Samples from true function
trPRS=0;                                                 % Ordinary Kriging
[XX,YY]=meshgrid(-5:dN:10, 0:dN:15);              % NxN grid of predictions
xp=[reshape(XX,N^2,1) reshape(YY,N^2,1)];               % Prediction points
[yp, ysig, nb]=Kriging(xs,ys,xp,trPRS,0);         % Build Kriging surrogate
YP=reshape(yp,N,N);                                   % Predictions at grid
YT=reshape(BraninHu(xp),N,N);                       % True function at grid
figure(1); scatter(xs(:,1),xs(:,2),'filled');       % Plot sample locations
axis([-5 10 0 15]); xlabel('x_1');ylabel('x_2');
figure(2); hold on; surf(XX,YY,YP);                   % Plot of predictions
plot3(xs(:,1),xs(:,2),ys,'ro');                           % Plot of samples
axis([-5 10 0 15]);xlabel('x_1');ylabel('x_2');zlabel('y');
figure(3); hold on; surf(XX,YY,YT)                  % PLot of true function
plot3(xs(:,1),xs(:,2),ys,'ro');                           % Plot of samples
axis([-5 10 0 15]);xlabel('x_1');ylabel('x_2');zlabel('y');
eRMS=sqrt(sum(sum((YT-YP).^2))/N^2)          % RMS error at NxN grid points
eav=sum(sum(abs(YT-YP)))/N^2             % Average error at NxN grid points
emax=max(max(abs(YT-YP)))                % Maximum error at NxN grid points
%-----------------------------------------------------------------------------
%
% Figure 7.26
BraninHu = @(x) (x(:,2)-1.275*(x(:,1)/pi).^2+5*(x(:,1)/pi)-6).^2 ...
                 +10*(1-1/(8*pi))*cos(x(:,1))+10;
rng('default'); ny=20;
xs=lhsdesign(ny,2,'iteration',5000);                       % 20 LHS samples
xs=xs.*[15, 15]+[-5, 0];                    %xs=[xs;-5 0;10 0;-5 15;10 15];
ys=BraninHu(xs);                               % Samples from true function
trPRS=0;                                                 % Ordinary Kriging
s=(0:0.01:1)'; N=size(s,1);          % Parameter along a diagonal direction
xp=[-5+15*s 15*s];                                  % xp=[-5+15*s 15-15*s];
[yp,ysig,nb]=Kriging(xs,ys,xp,trPRS,0);           % Build Kriging surrogate
yt=BraninHu(xp);                             % True function for comparison
cilb = yp + tinv(0.025,ny-nb)*ysig;  % Lower-bounds of confidence intervals
ciub = yp + tinv(0.975,ny-nb)*ysig;  % Upper-bounds of confidence intervals
figure; hold on;                                         % Plotting results
xShaded=[s;sort(s,'descend')]; yShaded=[cilb;ciub(end:-1:1)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(s,yp,'r-',s,yt,'k-','Markersize',10,'Linewidth',1.5)
xlabel('s');ylabel('y(s)');
legend('95% PI','Kriging prediction','True function');
eRMS=sqrt(sum((yt-yp).^2)/N);                % RMS error at NxN grid points
eav=sum(abs(yt-yp))/N;                   % Average error at NxN grid points
emax=max(abs(yt-yp));                    % Maximum error at NxN grid points
fprintf('RMS error:%8.3g\nAve error:%8.3g\nMax error:%8.3g\n',eRMS,eav,emax);
%-----------------------------------------------------------------------------
%
% Example 7.15
myfun = @(x) sin(x);                                     % In-line function
rng('default'); ny=10; nstd=0.1;            % No. of samples, and noise std
xs=linspace(0,2*pi,ny)';                                          % Samples
ys=myfun(xs);                                  % Samples from true function
noise=lhsnorm(0,nstd^2,ny); noise=nstd*(noise-mean(noise))/std(noise);
fprintf('Sample noiseSTD:%8.3g\n',std(noise));
ys=ys+noise;
trPRS=0;                                                 % Ordinary Kriging
xp=(0:0.01:2*pi)';                                      % Prediction points
[yp, ysig, nb]=Kriging(xs,ys,xp,trPRS,1);         % Build Kriging surrogate
cilb = yp + tinv(0.025,ny-nb)*ysig;  % Lower-bounds of confidence intervals
ciub = yp + tinv(0.975,ny-nb)*ysig;  % Upper-bounds of confidence intervals
ytrue=myfun(xp);                 % True function values at predictio points
figure; hold on;                                         % Plotting results
xShaded=[xp;sort(xp,'descend')]; yShaded=[cilb;ciub(end:-1:1)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(xp,yp,'r-',xs,ys,'k+',xp,ytrue,'k-','Markersize',10,'Linewidth',1.5)
xlabel('x');ylabel('y(x)');
legend('95% PI','Kriging prediction','Samples','True function');
%-----------------------------------------------------------------------------