function [CImin, CImax, CIx]=Confidence_bands_new_observation(x,y,yf,p)
%Finds values for the confidence band of a new observation
%
%Usage:
%[CImin, CImax, CIx]=Confidence_bands_new_observation(x,y,yf,p)
%
%Input:
%x   [n]   x values
%y   [n]   y values, usaually measured
%yf  [n]   fitted values, e.g. for model fitted trough x vs y relationships
%p         probability for exceeding confidence interval, e.g. 0.05 for 95% of obs 
%          in confidence interval
%
%Output
%CImin  [n] sorted y values for lower confidence band
%CImax  [n] sorted y values for upper confidence band 
%CIx    [n] sorted x values for confidence bands
%
%Tom Schut 2007

if nargin<4 p=0.9; end

a=find(~isnan(mean(x,2)) & ~isnan(y));
b=find(isnan(mean(x,2)) | isnan(y));
x=x(a,:); y=y(a); yf=yf(a);

nd=size(x,1);
[t]=tinv(1-p/2,nd);


ERR=(y-yf);
SE=ERR.^2;
SSX=(x-SE).^2;
SPE=(x-yf).^2;
MSE=sum(SE)./nd;
S=sqrt(MSE.*(1/nd + 1/1 + SSX./sum(SSX)));

[CIx ix]=sort(x);
CImin=yf(ix)-t.*S;
CImax=yf(ix)+t.*S;
