clear all, close all, clc


p = 10;                              % number of terms in the model
x = linspace(0,4,100)';              % number of observations
f = 3*x + 2*x.^2 + 5*randn(100,1);   % observations

A = zeros(100,p);
for j=1:p
    A(:,j) = x.^(j-1);
end


xL2 = pinv(A)*f;

[XL1 FitInfo] = lasso(A,f,'CV',10);
%lassoPlot(XL1,FitInfo,'PlotType','CV')
%lassoPlot(XL1,FitInfo,'PlotType','Lambda')

xL1 = XL1(:,FitInfo.Index1SE);
xL1DeBiased = pinv(A(:,abs(xL1)>0))*f;

%%%%

x1 = linspace(6,7,50);         % future data
f1 = 3*x1 + 2*x1.^2 + 5*randn(50,1);

% least-squares error
E2 = 0;
for k=1:50
    funcVal2 = 0;
    for j=1:p
        funcVal2 = funcVal2 + xL2(j)*x1(k)^(j-1);
    end
    E2 = E2 + abs( funcVal2 - f1(k) )^2;
end

% LASSO error
E1 = 0;
ind = find(abs(xL1)>0);
M = size(ind,1);
for k=1:50
    funcVal1 = 0;
    for j=1:M
        funcVal1 = funcVal1 + xL1DeBiased(j)*x1(k)^(ind(j)-1);
    end
    E1 = E1 + abs( funcVal1 - f1(k) )^2;
end


