close all; clear; clc;

%%
%rng(41);
seed = round(10000*rand);
rng(4089);

x = linspace(-1,1,100);
y = 3*x - 1 + 0.3*randn(size(x));
y(randi(50,1,10)) = 5*rand(1,10);
m = @(p) p(1)*x + p(2);

loss_L2 = @(p) mean( (m(p)-y).^2 ); % L2 norm
loss_L1 = @(p) mean( abs(m(p)-y) ); % L1 norm
loss_Linf = @(p) max( abs(m(p)-y) ); % Linf norm
% MLE with normal and extreme value noise models
loss_MLEnorm = @(p) -sum(log( normpdf(m(p(1:2))-y,0,p(3)) ));
loss_MLEext = @(p) -sum(log( pdf('Extreme Value',m(p(1:2))-y,0,p(3)) ));

pf_L2 = fminunc(loss_L2,[0 0]);
pf_L1 = fminunc(loss_L1,[0 0]);
pf_Linf = fminunc(loss_Linf,[0 0]);

V = [x(:) ones(size(x(:)))];
pf_linalg = V\y(:);

pf_MLEnorm = fminunc(loss_MLEnorm,[0 0 1]);
pf_MLEext = fminunc(loss_MLEext,[0 0 1]);

T = speye(numel(y),numel(y));

%% Plot
cm = lines;
figure(1); clf;
a = axes; hold on;
plot(0,0,'k+')
plot(x,y,'k.');
%plot(x,m([3,-1]),'k-');
plot(x,polyval(pf_L2,x));
plot(x,polyval(pf_L1,x));
plot(x,polyval(pf_Linf,x));
%plot(x,polyval(pf4(1:2),x),'--');
plot(x,polyval(pf_MLEext(1:2),x),'--');
plot(x,polyval(pf_linalg,x),'--');
xlim([-1 1]); ylim([-5 5]);
axis square; grid minor; box on;
xlabel('x'); ylabel('y');
a.XTickLabel = ''; a.YTickLabel = '';