%
%                   Example 13.4
%
%  Objective: Compare 2-sigma confidence intervals with frequentist prediction intervals
%
%  Notation: P_data, psi_data - Data in the design matrix X used to compute the mean
%                               and confidence intervals
%            P_test, psi_test - New observations used in X_test to compute prediction intervals
%            V - Covariance matrix for parameters
%            V_obs - Covariance matrix for observation errors
%

  clear all
  %close all

%
% Load the data on [0,0.6] and specify the parameter mean, covariance, 
% and observation error variance
%

  load Helmholtz2_new.txt
  P_data = Helmholtz2_new(:,1);
  psi_data = Helmholtz2_new(:,2);

  Pmin = 0;
  Pmax = 0.6;
  dP = 0.01;
  P_datat =  Pmin:dP:Pmax;
  P_data = P_datat';
  n = length(P_data);
  p = 3;
  X = [P_data.^2 P_data.^4 P_data.^6];
  H = inv(X'*X);

  sigma = 3.6;
  error = sigma*randn(size(P_data));

%%%% Construct synthetic data and save as Helmholtz2_new.txt.  Comment out the next
%    4 lines if using the data in Helmholtz2_new.txt. 

%  alpha_1 = -382.7;
%  alpha_11 = 760.3;
%  alpha_111 = 155.5;

%  alpha_1 = -482.7;
%  alpha_11 = 860.3;
%  alpha_111 = 155.5;

%  psi_data = alpha_1*P_data.^2 + alpha_11*P_data.^4 + alpha_111*P_data.^6 + error;

  q_bar = H*(X')*psi_data;
  alpha_1 = q_bar(1);
  alpha_11 = q_bar(2);
  alpha_111 = q_bar(3);

  psi = alpha_1*P_data.^2 + alpha_11*P_data.^4 + alpha_111*P_data.^6;
  R = psi_data - psi;
  sigma2 = (1/(n-p))*R'*R;

  q_bar = H*(X')*psi_data;
  alpha_1 = q_bar(1);
  alpha_11 = q_bar(2);
  alpha_111 = q_bar(3);

  P = P_data;
  psi_alpha_1 = P.^2;
  psi_alpha_11 = P.^4;
  psi_alpha_111 = P.^6;
  sensmat = [psi_alpha_1 psi_alpha_11 psi_alpha_111];
  V = sigma2*inv(sensmat'*sensmat);
  q_bar = H*(X')*psi_data;

  udata = psi_data;
  data = [P udata];
  %save('Helmholtz2_new.txt','data','-ascii') 

  format short e
  q_bar
  2*sqrt(V(3,3))

%%%%%%
  
  var = sigma^2;
  V_sig = (sigma^2)*eye(n,n);

%
% Compute the response mean, covariance and a 2-standard deviation bound
%

  mean = X*q_bar;
  covY = X*V*X' + V_sig;
  sd_bd = 2*sqrt(diag(covY));

%
% Construct the test values P_test and store in X_test. Compute the alpha/2 prediction
% intervals and plot the intervals for P in [0,0.6] and [0,0.8].
%

  P_test = Pmin:dP:Pmax;
  psi_test = alpha_1*P_test.^2 + alpha_11*P_test.^4 + alpha_111*P_test.^6;
  N = length(P_test);
  X_test = [P_test'.^2 P_test'.^4 P_test'.^6];

  alpha = 0.05;
  tcrit2 = tinv(1-alpha/2,n-p);
  for j=1:N
    int_upper(j) = psi_test(j) + tcrit2*sigma*sqrt(1 + X_test(j,:)*inv(X'*X)*(X_test(j,:)'));
    int_lower(j) = psi_test(j) - tcrit2*sigma*sqrt(1 + X_test(j,:)*inv(X'*X)*(X_test(j,:)'));
  end

  figure(1)
  plot(P_test,psi_test,'-k','linewidth',3)
  hold on
  plot(P_data,psi_data,'xr','linewidth',3)
  plot(P_data,mean+sd_bd,'-.k','linewidth',3)
  plot(P_test,int_upper,'--b','linewidth',3)
  plot(P_data,mean-sd_bd,'-.k','linewidth',3)
  plot(P_test,int_lower,'--b','linewidth',3)
  hold off
  set(gca,'FontSize',[20])
  legend('Mean','Data','2\sigma Interval','Prediction Interval','Location','NorthEast')
  xlabel('Polarization P')
  ylabel('Helmholtz Energy \psi')
  %print -depsc Polarization_pred_ints_a

  dP = 0.01;
  Pmin = 0;
  Pmax = 0.8;
  P_test = Pmin:dP:Pmax;
  psi_test = alpha_1*P_test.^2 + alpha_11*P_test.^4 + alpha_111*P_test.^6;
  N = length(P_test);
  X_test = [P_test'.^2 P_test'.^4 P_test'.^6];

  alpha = 0.05;
  tcrit2 = tinv(1-alpha/2,n-p);
  for j=1:N
    int_upper(j) = psi_test(j) + tcrit2*sigma*sqrt(1 + X_test(j,:)*inv(X'*X)*(X_test(j,:)'));
    int_lower(j) = psi_test(j) - tcrit2*sigma*sqrt(1 + X_test(j,:)*inv(X'*X)*(X_test(j,:)'));
  end

  figure(2)
  plot(P_test,psi_test,'-k','linewidth',3)
  hold on
  plot(P_data,psi_data,'xr','linewidth',3)
  plot(P_test,int_upper,'--b','linewidth',3)
  plot(P_test,int_lower,'--b','linewidth',3)
  hold off
  set(gca,'FontSize',[20])
  legend('Mean','Data','Prediction Interval','Location','NorthWest')
  xlabel('Polarization P')
  ylabel('Helmholtz Energy \psi')
  %print -depsc Polarization_pred_ints_b