% Demo prepared by Romain Thibaux <thibaux@eecs.berkeley.edu>
% for the Regression lecture of the Practical Machine Learning
% class (CS294 - Fall 2006).

% Script for lecture 2 of the Practical ML class
P = dlmread('mote_locs.txt');
P = P(:,2:end);
N = size(P,1);

X1 = dlmread('lab_temperature_1.data');
scatter3(P(:,1),P(:,2),X1(:,1),'r','filled');
axis vis3d;
title('Temperature at noon', 'FontSize', 18);
% zlabel('Temperature', 'FontSize', 18);
% xlabel('x', 'FontSize', 18);
% ylabel('y', 'FontSize', 18);
pause;
I = [1:4, 6:N];  % 5 is faulty
scatter3(P(I,1),P(I,2),X1(I,1),'r','filled');
axis vis3d;
title('Temperature at noon (sensor 5 removed)', 'FontSize', 18);
pause;

PP = [ones(N,1), P];
%w = PP(I,:) \ X1(I,1);
w = [PP'*PP]\[PP'*X1(:,1)];
x = linspace(min(P(:,1)),max(P(:,1)),10)';
y = linspace(min(P(:,2)),max(P(:,2)),10)';
X = ones(size(y))*x';
Y = y*ones(size(x'));
Z = w(1) + w(2)*X + w(3)*Y;
hold on;
surf(x, y, Z);
title('Temperature at noon with linear fit', 'FontSize', 18);
hold off;
pause;

X2 = dlmread('lab_temperature_3.data');
scatter3(P(I,1),P(I,2),X2(I,1),'r','filled');
axis vis3d;
title('Temperature at midnight', 'FontSize', 18);
pause;

PP = [ones(N,1), P, P.^2];
PP = PP(I,:);
w = [PP'*PP]\[PP'*X2(I,1)];
Z = w(1) + w(2)*X + w(3)*Y + w(4)*X.^2 + w(5)*Y.^2;
hold on;
surf(x, y, Z);
title('Temperature at midnight with quadratic fit', 'FontSize', 18);
zlabel('Temperature', 'FontSize', 18);
xlabel('x', 'FontSize', 18);
ylabel('y', 'FontSize', 18);
axis manual;
hold off;
pause;

scatter3(P(:,1),P(:,2),X2(:,1),'r','filled');
hold on;
surf(x, y, Z);
axis vis3d;
title('Temperature at midnight with quadratic fit', 'FontSize', 18);
hold off;
pause;

T = view;
scatter3(P(:,1),P(:,2),X2(:,1),'r','filled');
PP = [ones(N,1), P, P.^2];
w = [PP'*PP]\[PP'*X2(:,1)];
Z = w(1) + w(2)*X + w(3)*Y + w(4)*X.^2 + w(5)*Y.^2;
hold on;
surf(x, y, Z);
title('Temperature at midnight with quadratic fit (including sensor 5)', 'FontSize', 18);
axis vis3d;
view(T);
hold off;
pause;

randn('state', 654321);
w = [-1.5; 1/9];
y = w(1)*(1:20)+w(2)*(1:20).^2;
z = y+2*randn(1,20);
plot(1:20,y,'k');
hold on;
scatter(1:20,z,'r','filled');
hold off;
title('Simple quadratic fit', 'FontSize', 18);
pause;

scatter3(1:20,(1:20).^2,z,'b','filled');
hold on;
m = min(z);
M = max(z);
line(1:20, (1:20).^2, m*ones(1,20), 'LineWidth', 2);
line(1:20, (1:20).^2, M*ones(1,20), 'LineWidth', 2);
line([1;1],[1;1],[m,M], 'LineWidth', 2);
line([20;20],[400;400],[m,M], 'LineWidth', 2);
K = [20, 400];
k = [-400, 20];
k = k/norm(k);
k = 2*max(abs(k*[1:20;(1:20).^2]))*k;
title('Geometric interpretation', 'FontSize', 18);
xlabel('x', 'FontSize', 18);
ylabel('x^2', 'FontSize', 18);
X = [k(1); k(1)+K(1); -k(1)+K(1); -k(1); k(1)];
Y = [k(2); k(2)+K(2); -k(2)+K(2); -k(2); k(2)];
Z = [k*w; (k+K)*w; (-k+K)*w; -k*w; k*w];
axis([min(X),max(X),min(Y),max(Y),min(Z),max(Z)]);
axis vis3d;
pause;

line(X, Y, Z, 'Color', 'r', 'LineWidth', 2);
plot3(1:20,(1:20).^2,y,'r', 'LineWidth', 2);
hold off;
pause;

V = 1:0.1:20;
X = zeros(20,0);
XX = zeros(length(V),0);
for i=0:19
  X = [X, ((-9:10)/10)'.^i];
  XX = [XX, ((-9:0.1:10)/10)'.^i];
  w = (X'*X)\(X'*z');
  Y = XX*w;
  scatter(1:20,z,'b','filled');
  hold on;
  plot(V,Y);
  hold off;
  title(sprintf('Runge''s phenomenon (degree %d)', i), 'FontSize', 18);
  pause;
end

eps = 0.1;
V = 1:0.1:20;
X = zeros(20,0);
XX = zeros(length(V),0);
for i=0:19
  X = [X, ((-9:10)/10)'.^i];
  XX = [XX, ((-9:0.1:10)/10)'.^i];
  w = (X'*X + eps*eye(i+1))\(X'*z');
  Y = XX*w;
  scatter(1:20,z,'b','filled');
  hold on;
  plot(V,Y);
  hold off;
  title(sprintf('Effect of regularization (degree %d)', i), 'FontSize', 18);
  pause;
end

for sigma = [10,4,2,1,0.5,0.2,0.1]
  Kernel = exp(-0.5*(-30:0.1:30).^2/(sigma^2));
  L = (length(Kernel)+1)/2;  % midpoint
  for i=1:length(V)
    k = Kernel(L+(0:19)*10-i);
    Y(i) = k*z'/sum(k);
  end
  scatter(1:20,z,'b','filled');
  hold on;
  plot(V,Y);
  title(sprintf('Kernel regression (sigma=%d)', sigma), 'FontSize', 18);
  hold off;
  pause;
end