%
%                      Example18_6.m 
%

  clear all
  %close all

%
% Create training points and responses along with the test points and their responses
%

  a = 0.05;   
  b = 0.05;  

  xs = lhsdesign(10,1);
  Y = @(x) ((x-1/4).^2 + 1/2 + a*(x+b).*sin(20*pi*x))';
  ys = Y(xs);
  xp = linspace(0,1,100)';
  yp = Y(xp);

%
% Input initial values for hyperparameters 
%   Book:   ell_k,  sigma,  sigma_0
%   MATLAB: sigmaL, sigmaF, sigma 
%

  kparams = [3.5, 10];          % [sigmaL, sigmaF]

%
% Estimate the hyperparameters using an initial estimate of sigma = eps and compute the expected
% prediction at xp along with the variance
%

  gprMdl = fitrgp(xs,ys,'KernelFunction','squaredexponential','KernelParameters',kparams,'Sigma',eps);
  [pred,~,yint] = predict(gprMdl,xp);

%
%  Extract the optimized values of sigma, sigmaL, and sigmaF
%  Extract the covariance function for plotting
%

  sigmaL = gprMdl.KernelInformation.KernelParameters(1); 
  sigmaF = gprMdl.KernelInformation.KernelParameters(2); 
  sigma  = gprMdl.Sigma; 
  kfcn = gprMdl.Impl.Kernel.makeKernelAsFunctionOfXNXM(gprMdl.Impl.ThetaHat);
  K = kfcn(xp(1),xp(1:end));

%
% Plot the covariance function and predictions and standard deviation intervals.
%

  figure(1)
  plot(xp,K,'b-','linewidth',3)
  %axis([0 1.5 0 0.2])
  set(gca,'Fontsize',[22]);
  xlabel('Parameter q') 
  ylabel('Covariance Function c')


  figure(2)
  f = [yint(:,2); flipud(yint(:,1))];
  h(1) = fill([xp; flipud(xp)], f, [7 7 7]/8)
  set(get(get(h(1),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
  hold on
  h(2) = plot(xs,ys,'ro','linewidth',5,'DisplayName','Data');
  h(3) = plot(xp,pred,'b-','linewidth',3,'DisplayName','Predictive Mean');
  hold off
  legend('Location','NorthWest')
  set(gca,'Fontsize',[22]);
  xlabel('Parameter q')
  ylabel('Response')