%
%        Example13_9.m
%


  clear all

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

  Pmin = 0;
  Pmax = 0.6;
  Pmax2 = 0.8;
  dP = 0.01;

%
% This next bit should be commented out later.
%
%  P_data = [0:0.01:0.6]';
  sigma = 5.6;
  error = sigma*randn(size(P_data));
  alpha_1 = -382.19;
  alpha_11 = 713.46;
  alpha_111 = 150.50;
%  psi_data = alpha_1*P_data.^2 + alpha_11*P_data.^4 + alpha_111*P_data.^6 + error;
%

  n = length(P_data);
  p = 3;
  X = [P_data.^2 P_data.^4 P_data.^6];
  H = inv(X'*X);

  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;
  sigma1 = sqrt(sigma2);
  alpha = 0.05;
  tcrit2 = tinv(1-alpha/2,n-p);

  V = sigma2*H;

  for j=1:81
    p0 = (j-1)*dP;
    P0(j) = p0;
    pp0 = [p0^2; p0^4; p0^6];
    factor1 = sqrt(pp0'*H*pp0);
    factor2 = sqrt(1 + pp0'*H*pp0);
    psi_mean = alpha_1*p0^2 + alpha_11*p0^4 + alpha_111*p0^6;
    Psi_mean(j) = psi_mean;
    conf1(j) = psi_mean - tcrit2*sigma1*factor1;
    conf2(j) = psi_mean + tcrit2*sigma1*factor1;
    pred1(j) = psi_mean - tcrit2*sigma1*factor2;
    pred2(j) = psi_mean + tcrit2*sigma1*factor2;
  end

  figure(1)
  box on
  dimc = [0.90, 0.90, 0.90];
  fillyy(P0,pred1,pred2,dimc)
  hold on
  dimc = [0.70, 0.70, 0.70];
  fillyy(P0,conf1,conf2,dimc)
  plot(P0,Psi_mean,'-k','linewidth',2)
  plot(P_data,psi_data,'xk','linewidth',3)
%  axis([0 0.8 -100 150])
  hold off
  set(gca,'Fontsize',[20]);
  xlabel('Polarization P')
  ylabel('Helmholtz Energy \psi')
%  print -depsc intervals_frequentist

%
% Construct DRAM parameters and settings
%

  clear data model options

  data.xdata = P_data;
  data.ydata = psi_data;
  tcov = V;
  tmin = [alpha_1; alpha_11; alpha_111];

  params = {
    {'q1',tmin(1), -800, -100}
    {'q2',tmin(2), 0, 2000}
    {'q3',tmin(3), 0, 2000}
  };
  model.ssfun = @landauss;
  model.sigma2 = sigma2;
  options.qcov = tcov;
  options.nsimu = 20000;
  options.updatesigma = 1;

%
% Run DRAM to construct the chains (stored in chain) and measurement
% variance, which is stored in s2chain.
%

  [results,chain,s2chain] = mcmcrun(model,data,params,options);

%
%
  figure(3)
  mcmcplot(chain,[],results,'chainpanel');

  figure(4)
  mcmcplot(chain,[],results,'pairs');

%
% Compute credible and prediction intervals
%
 
  xvals = linspace(0,Pmax2)';
  nsample = 5000;
  out = mcmcpred(results,chain,s2chain,xvals,@energy_eval,nsample)


  alpha_1_val = chain(10:10:end,1);
  alpha_11_val = chain(10:10:end,2);
  alpha_111_val = chain(10:10:end,3);
  error_end = sigma1*randn(size(alpha_1_val));

  psi_end = alpha_1_val*Pmax2^2 + alpha_11_val*Pmax2^4 + alpha_111_val*Pmax2.^6 + error_end;
  [bandwidth_psi_end,density_psi_end,psi_end_mesh,cdf_psi_end]=kde(psi_end);
  nbins = 20;

  lower_lim = prctile(psi_end,2.5);
  upper_lim = prctile(psi_end,97.5);

  pvals = [0.025,0.975];
  [n,m] = size(psi_end); if n==1; n=m;end
  limits = interp1(sort(psi_end),(n-1)*pvals+1);


  figure(4)
  mcmcpredplot2(out);
%  mcmcpredplot(out);
  hold on
  box on
  plot(P_data,psi_data,'xk','linewidth',3)
  plot([0.77 0.8],[lower_lim lower_lim],'r-','linewidth',3)
  plot([0.77 0.8],[upper_lim upper_lim],'r-','linewidth',3)
%  axis([0 Pmax2 -100 150])
  hold off
  set(gca,'Fontsize',[20]);
  xlabel('Polarization P')
  ylabel('Helmholtz Energy \psi')
  %print -depsc intervals_Bayesian


  figure(5)
%  histnorm(psi_end,nbins);
  plot(psi_end_mesh,density_psi_end,'b-','linewidth',3)
  y1 = get(gca,'ylim');
  hold on
  plot([lower_lim lower_lim],y1,'r--','linewidth',3)
  plot([upper_lim upper_lim],y1,'r--','linewidth',3)
  hold off
  set(gca,'Fontsize',[20]);
  xlabel('Helmholtz Energy \psi(0.8)')
  ylabel('Posterior Predictive Density')
  %print -depsc predictive_density