function h=mcmcpredplot(out,data,adddata)
%MCMCPREDPLOT - predictive plot for mcmc results
% Creates predictive figures for each batch in the data set using
% mcmc chain. Needs input from the function mcmcpred.
% Example:
%  out=mcmcpred(results,chain,s2chain,data,modelfun);
%  mcmcpredplot(out)
%
% If s2chain has been given to mcmcpred, then the plot shows 95%
% probability limits for new observations and for model parameter
% uncertainty. If s2chain is not used then the plot contains 50%,
% 90%, 95%, and 99% predictive probability limits due parameter uncertainty.

% $Revision: 1.4 $  $Date: 2007/08/22 16:10:58 $

if nargin < 2
  data = out.data;
end

if nargin < 3
  adddata = 0; % add data.ydata datapoints to the figures
end

nbatch = length(out.predlims);

if ~iscell(data)
  d=data; data=[]; data{1}=d; clear d
end

np = size(out.predlims{1}{1},1);
nn = (np+1)/2; % median
np = nn-1;

hh = zeros(nbatch,1);

for i=1:nbatch
  if nbatch > 1; hh(i) = figure;else hh(1)=gcf; end; % create new figures
  plimi = out.predlims{i};
  ny = size(plimi,2);
  
  datai = data{i};
  
  if isnumeric(datai)
    time = datai(:,1); % time is the first columd of data
  elseif isfield(datai,'ydata')
    time = datai.ydata(:,1); % first column of ydata
  elseif isfield(datai,'xdata')
    time = datai.xdata(:,1); % first column of xdata
  else
    error('dont know the x axis of the plots')
  end
    
  for j=1:ny
    dimc = [0.9 0.9 0.9]; % dimmest (lightest) color
    dimc = [0.75 0.75 0.75];
    %subplot(ny,1,j)
    figure
    set(gca,'Fontsize',[22]);
    if ~isempty(out.obslims)
      line1 = out.obslims{i}{j}(1,:);
      line2 = out.obslims{i}{j}(3,:);
      %fillyy(time,out.obslims{i}{j}(1,:),out.obslims{i}{j}(3,:),dimc);
      fillyy(time,line1,line2,dimc);
      hold on
      %dimc = [0.8,0.8,0.8];
      dimc = [0.5, 0.5, 0.5];
    end
    fillyy(time,plimi{j}(1,:),plimi{j}(2*nn-1,:),dimc);
    hold on
    for k=2:(nn-1)
      fillyy(time,plimi{j}(k,:),plimi{j}(2*nn-k,:),dimc.*0.9.^(k-1));
    end
    plot(time,plimi{j}(nn,:),'-k','linewidth',2);
    if adddata
      plot(datai.ydata(:,1), datai.ydata(:,j+1),'sk');
      xlabel('Time (s)')
    end
    hold off
    if nbatch > 1
%      title(sprintf('Data set %d, y[%d]',i,j));
    elseif ny > 1
%      title(sprintf('y[%d]',j));
    end
  end
  
end

if nargout > 0
  h=hh;
end