clear all
close all
fNum = 1;

%% Demo of Least Squares
nInstances = 500;
X = randn(nInstances,1);
y = X*(1/2) +randn(nInstances,1)/4; % y is bilinear in X(i,:) and w, plus random Gaussian noise

% Least Squares
figure(fNum);fNum=fNum+1;clf; hold on;
plot(X,y,'.')
wLS = (X'*X)\(X'*y);
domain = xlim;
domain = domain(1):.01:domain(2);
for i = 1:length(domain)
    f(i) = wLS*domain(i);
end
plot(domain,f,'r');

%% Demo of Basis Expansion

% Generate polynomial data
y = (-X.^2)*(1/2) + 1 + randn(nInstances,1)/4;

figure(fNum);fNum=fNum+1;clf; hold on;
plot(X,y,'.');

% Least Squares
wLS = (X'*X)\(X'*y);
domain = xlim;
domain = domain(1):.01:domain(2);
for i = 1:length(domain)
    f(i) = wLS*domain(i);
end
plot(domain,f,'r');


figure(fNum);fNum=fNum+1;clf; hold on;
plot(X,y,'.');

% Least Squares on Expanded Basis
Xbasis = [ones(nInstances,1) X X.^2];
wLSbasis = (Xbasis'*Xbasis)\(Xbasis'*y);
for i = 1:length(domain)
    fbasis(i) = wLSbasis'*[1;domain(i);domain(i)^2];
end
plot(domain,fbasis,'g');

%% Demo of Overfitting and Regularization
nInstances = 9;
X = randn(nInstances,1);
y = (X.^4)*(1) + (X.^3)*(-3) + (X.^2)*3 + X*-3 + 1 + randn(nInstances,1)/4;
y = y-mean(y);

figure(fNum);fNum=fNum+1;clf; hold on;
plot(X,y,'.');

% Least Squares
wLS = (X'*X)\(X'*y);
domain = xlim;
domain = domain(1):.01:domain(2);
for i = 1:length(domain)
    f1(i) = wLS*domain(i);
end
plot(domain,f1,'r');

figure(fNum);fNum=fNum+1;clf; hold on;
plot(X,y,'.');

% Degree 2 Polynomial
Xbasis = [ones(nInstances,1) X X.^2];
wLSbasis = (Xbasis'*Xbasis)\(Xbasis'*y);
for i = 1:length(domain)
    f2(i) = wLSbasis'*[1;domain(i);domain(i)^2];
end
plot(domain,f2,'g');
yl = ylim;

figure(fNum);fNum=fNum+1;clf; hold on;
plot(X,y,'.');

% Degree 4 Polynomial
Xbasis = [ones(nInstances,1) X X.^2 X.^3 X.^4];
wLSbasis = (Xbasis'*Xbasis)\(Xbasis'*y);
for i = 1:length(domain)
    f4(i) = wLSbasis'*[1;domain(i);domain(i)^2;domain(i)^3;domain(i)^4];
end
plot(domain,f4,'g');
ylim(yl);

figure(fNum);fNum=fNum+1;clf; hold on;
plot(X,y,'.');

% Degree 8 Polynomial
Xbasis = [ones(nInstances,1) X X.^2 X.^3 X.^4 X.^5 X.^6 X.^7 X.^8];
wLSbasis = (Xbasis'*Xbasis)\(Xbasis'*y);
for i = 1:length(domain)
    f8(i) = wLSbasis'*[1;domain(i);domain(i)^2;domain(i)^3;domain(i)^4;domain(i)^5;domain(i)^6;domain(i)^7;domain(i)^8];
end
plot(domain,f8,'g');
ylim(yl);

figure(fNum);fNum=fNum+1;clf; hold on;
plot(X,y,'.');

% Regularized Degree 8 Polynomial
Xbasis = [ones(nInstances,1) X X.^2 X.^3 X.^4 X.^5 X.^6 X.^7 X.^8];
regularizer = 10*eye(9);
regularizer(1) = 0;
wLSbasis = (Xbasis'*Xbasis + regularizer)\(Xbasis'*y);
for i = 1:length(domain)
    f8R(i) = wLSbasis'*[1;domain(i);domain(i)^2;domain(i)^3;domain(i)^4;domain(i)^5;domain(i)^6;domain(i)^7;domain(i)^8];
end
plot(domain,f8R,'g');
ylim(yl);

%% Linear vs. Sigmoid Output
domain = -5:.1:5;
for i = 1:length(domain)
    xLin(i) = domain(i);
    xSig(i) = 1/(1+exp(-domain(i)));
end
figure(fNum);fNum=fNum+1;clf; hold on;
plot(domain,xLin);
yl = ylim;
figure(fNum);fNum=fNum+1;clf; hold on;
plot(domain,xSig);
ylim(yl);

%% Loss Functions
domain = -2:.1:2;
epsilon = .5;
for i = 1:length(domain)
    errLS(i) = domain(i)^2;
    errL1(i) = abs(domain(i));
    errSV(i) = max(abs(domain(i))-epsilon,0);
end
figure(fNum);fNum=fNum+1;clf; hold on;
plot(domain,zeros(size(domain)),'k-');
plot(domain,errLS);
yl = ylim;
yl(1) = -1;
ylim(yl);
figure(fNum);fNum=fNum+1;clf; hold on;
plot(domain,zeros(size(domain)),'k-');
plot(domain,errL1);
ylim(yl);
figure(fNum);fNum=fNum+1;clf; hold on;
plot(domain,zeros(size(domain)),'k-');
plot(domain,errSV);
ylim(yl);

%% Support Vector Regression
svrLoss = @(w,X,y,epsilon) sum(max(abs(X*w-y)-epsilon,0).^2);
options.numDiff = 1;
epsilon = 1;

nInstances = 200;
X = randn(nInstances,1);
y = (-X.^2)*(1/2) + 1 + randn(nInstances,1)/2;

figure(fNum);fNum=fNum+1;clf; hold on;
plot(X,y,'.');

Xbasis = [ones(nInstances,1) X X.^2];
wSVR = minFunc(svrLoss,zeros(size(Xbasis,2),1),options,Xbasis,y,epsilon);
domain = xlim;
domain = domain(1):.01:domain(2);
for i = 1:length(domain)
    fSVR(i) = wSVR'*[1;domain(i);domain(i)^2];
end
upper = fSVR+epsilon;
lower = fSVR-epsilon;
plot(domain,fSVR,'g');
plot(domain,upper,'c--');
plot(domain,lower,'c--');

% Find Support Vectors
SV = abs(Xbasis*wSVR - y) >= epsilon;
plot(X(SV==1),y(SV==1),'ro');

