%
%                        Example11_18.m  
%
%
% Construct prediction intervals for a quadratic problem and illustrate exptrapolation.
%
% Load the data and specify parameter values. The data was constructed using the commented
% commands.
%
  clear all

  load Y_data_pred

%  theta1 = 0.6;
%  theta2 = 1.2;
%  theta = [theta1 theta2];
%  sigma0 = 3.0;
%  x = 1:.1:3.5;
%  error = sigma0*randn(size(x));
%  Y = theta1 + theta2*x.^2 + error;
 
%
% Construct the data and compute parameter values based on this data. Compute the 95% tval. 
%

  x = 1:.1:3.5;
  X = [ones(size(x))' x'.^2];
  Finv = inv(X'*X);

  theta = Finv*(X')*Y';
  theta1 = theta(1,1);
  theta2 = theta(2,1);
  y = theta1 + theta2*x.^2;
  n = length(x);
  p = 2;
  R = Y - y;
  sigma2 = (1/(n-p))*R*R';
  sigma1 = sqrt(sigma2)
  alpha = 0.05;
  tval = tinv(1-alpha/2,n-p);   % tval = 2.0639;  95% Interval

  V = sigma2*Finv;

%
%  Construct confidence and prediction intervals on the interval [0,6]. 
%

  for j = 1:61
    x0 = (j-1)*0.1;
    X0(j) = x0;
    xx0 = [1; x0^2];
    factor1 = sqrt(xx0'*Finv*xx0);
    factor2 = sqrt(1 + xx0'*Finv*xx0); 
    Yh = theta1 + theta2*x0.^2;
    Yhvals(j) = Yh;
    conf1(j) = Yh - tval*sigma1*factor1;
    conf2(j) = Yh + tval*sigma1*factor1;
    pred1(j) = Yh - tval*sigma1*factor2;
    pred2(j) = Yh + tval*sigma1*factor2;
  end

%
% Plot the results.
% Figure 1: Different linetypes for confidence and prediction intervals
% Figure 2: Shaded regions for prediction and confidence intervals
% 

  figure(1)
  plot(X0,pred1,'--b','linewidth',3)
  hold on
  plot(X0,conf1,'-.k','linewidth',3)
  plot(X0,Yhvals,'-k','linewidth',3)
  plot(x,Y,'sr','linewidth',2)
  plot(X0,conf2,'-.k','linewidth',3)
  plot(X0,pred2,'--b','linewidth',3)
  hold off 
  set(gca,'Fontsize',[20]);
  legend('95% Prediction Inverval','95% Confidence Interval','Point Estimates','Data','Location','NorthWest')
  xlabel('x')
  ylabel('y')

  figure(2)
  box on
  dimc = [0.90, 0.90, 0.90];
  fillyy(X0,pred1,pred2,dimc)
  hold on
  dimc = [0.70, 0.70, 0.70];
  fillyy(X0,conf1,conf2,dimc)
  plot(X0,Yhvals,'-k','linewidth',2)
  plot(x,Y,'xb','linewidth',3)
  axis([0 6 -10 70])
%  xticks([0 1 2 3 4 5 6])
%  xticklabels({'0','1','2','3','4','5','6'})
  set(gca,'Fontsize',[22]);
  hold off
  xlabel('x')
  ylabel('y')