%
% problem 268
% example epage 250 
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');
load insulate; % befinsul

figure; 
h1=plot( befinsul(:,1), befinsul(:,2), 'rx' ); hold on; 
h2=plot( aftinsul(:,1), aftinsul(:,2), 'go' ); 
legend( [h1,h2], {'before','after'} ); 

% comment/uncomment the data set you want to run: 
% 
%x = befinsul(:,1); y = befinsul(:,2); tt='before insulation';
x = aftinsul(:,1); y = aftinsul(:,2); tt='after insulation';

fh=figure; ph1=plot( x, y, 'x' ); grid on; 
xlabel( 'average outside temperature' ); ylabel( 'weekly gas consumption per week per 1000' ); 
title( tt );

% cross validation with K=1
% (three polynomial models)
n = length(x); 
r1 = zeros(1,n); 
r2 = zeros(1,n); 
r3 = zeros(1,n); 
for i=1:n,
  xtest = x(i); 
  ytest = y(i); 
  xtrain = x; 
  ytrain = y; 
  xtrain(i)=[]; 
  ytrain(i)=[]; 
  [p1,s] = polyfit(xtrain,ytrain,1);  % first degree polynomial 
  [p2,s] = polyfit(xtrain,ytrain,2);  % second degree polynomial 
  [p3,s] = polyfit(xtrain,ytrain,3);  % third degree polynomial 
  r1(i) = (ytest-polyval(p1,xtest)).^2; 
  r2(i) = (ytest-polyval(p2,xtest)).^2; 
  r3(i) = (ytest-polyval(p3,xtest)).^2; 
end
% get the total prediction error for each model: 
pe1 = mean(r1); 
pe2 = mean(r2); 
pe3 = mean(r3); 

fprintf('%10d fold cross validation => estimated prediction error...\n',n); 
fprintf('(pe1,pe2,pe3) = (%10.5f,%10.5f,%10.5f)\n',pe1,pe2,pe3); 


% cross validation with K=4
% 
K = 4; 
r1 = zeros(1,n); 
r2 = zeros(1,n); 
r3 = zeros(1,n); 
r = floor(n/K);  % <- the number of samples per jackknife estimate ... rounded down 
for i=1:K,
  indsTest  = (r*(i-1)+1):(r*i);
  indsTrain = setdiff( 1:n, indsTest ); 
  xtest  = x( indsTest );
  ytest  = y( indsTest ); 
  xtrain = x( indsTrain ); 
  ytrain = y( indsTrain ); 
  [p1,s] = polyfit(xtrain,ytrain,1);  % first degree polynomial 
  [p2,s] = polyfit(xtrain,ytrain,2);  % second degree polynomial 
  [p3,s] = polyfit(xtrain,ytrain,3);  % third degree polynomial 
  r1(indsTest) = (ytest-polyval(p1,xtest)).^2; 
  r2(indsTest) = (ytest-polyval(p2,xtest)).^2; 
  r3(indsTest) = (ytest-polyval(p3,xtest)).^2; 
end
% remove the data points not used (because they don't fit into a multiple of K)
nf = r*K; 
r1 = r1(1:nf); r2 = r2(1:nf); r3 = r3(1:nf); 

% get the total prediction error for each model: 
pe1 = mean(r1); 
pe2 = mean(r2); 
pe3 = mean(r3); 

fprintf('%10d fold cross validation => estimated prediction error...\n',K); 
fprintf('(pe1,pe2,pe3) = (%10.5f,%10.5f,%10.5f)\n',pe1,pe2,pe3); 

% display the linear model
% 
figure(fh); hold on;
[p1,s] = polyfit(x,y,1); % fit to ALL the data 
%xplt=x.*(1+(1e-1)*randn(size(x))); ph2 = plot(xplt,polyval(p2,xplt),'-rx'); 
xplt=x; ph2 = plot(xplt,polyval(p1,xplt),'-rx'); 
legend( [ph1,ph2], {'raw data', 'first order poly'} ); 
saveas( gcf, 'prob_7_1_binsu_p1', 'epsc' ); 

pause; 

% display the quadratic model
% 
figure(fh); delete(ph2); hold on; 
[p2,s] = polyfit(x,y,2); % fit to ALL the data 
%xplt=x.*(1+(1e-1)*randn(size(x))); ph2 = plot(xplt,polyval(p2,xplt),'-rx'); 
xplt=x; ph2 = plot(xplt,polyval(p2,xplt),'-rx'); 
legend( [ph1,ph2], {'raw data', 'second order poly'} ); 
saveas( gcf, 'prob_7_1_binsu_p2', 'epsc' );