% This script demonstrates using a covariance matrix for to calculate a
% predictive interval.

%% Create a data set

% the model is a simple line

A=4.5;
B=0.5;

x=(1:10)';
y=A+B*x;

% add noise
sigma=5; %this is the sigma in normally distributed noise model
y_noise=y+sigma*randn(size(y));

%% Find A and B using Polyfit

n=1;
[p S]=polyfit(x,y_noise,n);
y_fit=polyval(p,x);

%Using our calcualted covariance
I=1/sigma^2*[sum(x.^2) sum(x); sum(x) length(x)];
V=inv(I);
y_1s=sqrt( V(1,1)*x_p.^2 + V(2,2)+2*V(2,1)*x_p);
y_1s_noCov=sqrt( V(1,1)*x_p.^2 + V(2,2)); %ignoring covariance (not correct)

%Using MATLAB Polyval predictive interval 
%Note MATLAB doc says it gives a 50% interval (not 1 sigma). 
x_p=(-10:50);
Rinv=inv(S.R)
normr=S.normr;
df=S.df;
V_pv=(Rinv*Rinv')*normr^2/df;
[y_p delta]=polyval(p,x_p,S);


%% Plot the results

figure
hold on
plot(x,y,'k','linewidth',2)
errorbar(x,y_noise,x*0+sigma,'r.','linewidth',2)
plot(x,y_fit,'g','linewidth',2)
plot(x_p,y_p+y_1s,'c','linewidth',2)
plot(x_p,y_p+y_1s_noCov,'m','linewidth',2)
plot(x_p,y_p-y_1s,'c','linewidth',2)
plot(x_p,y_p-y_1s_noCov,'m','linewidth',2)

legend('underlying model','noisy data','fit results','1\sigma interval',...
    'incorrect neglect of covariance')
xlabel('X ','fontsize',16)
ylabel('Y ','fontsize',16)
set(gca,'fontsize',14)


%% Test the predictive interval. 

figure
hold on
plot(x_p,y_p+y_1s,'c','linewidth',2)
plot(x_p,y_p-y_1s,'c','linewidth',2)

clear y_fittest
ytest=p(2)+x*p(1); %This assumes the estimated A,B are the true model
%In that case, the 1 sigma interval around this value will contain ~68% of future
%predictions.

for ii=1:100
y_noisetest=ytest+sigma*randn(size(ytest));
[p S]=polyfit(x,y_noisetest,n);
y_fittest(:,ii)=polyval(p,x_p);
plot(x_p,y_fittest(:,ii),'r','linewidth',0.5)
%pause(.1)
end

legend('1\sigma interval','1\sigma interval','fit results')
xlabel('X ','fontsize',16)
ylabel('Y ','fontsize',16)
set(gca,'fontsize',14)

%calculate some statistics
tst=40;
plot([x_p(40) x_p(40)],[ min(min(y_fittest)) max(max(y_fittest))],'k--')
frac=1/length(y_fittest(30,:))*sum(y_fittest(30,:)>(y_p(30)-y_1s(30))&...
    y_fittest(30,:)<(y_p(30)+y_1s(30)))









