% klagenfurt lectures collection
%
% MAXIMUM LIKELIHOOD ESTIMATION
%

%$ loading the data
load weightHeight.dat;    % the variable 
d_w = weightHeight(:,1);
d_h = weightHeight(:,2);
clear weightHeight;
Y = d_h;
sigma = 1000000;

%% visualisation init
xa = (42:3:120)';
outFormat = '%4d\t%4.2f\n';
clf; hold on; box on;
plot(d_w,d_h,'ro');
xlim([45 115]); ylim([135 195]);

%% MAXLIK using linear basis
X  = [d_w ones(size(d_h))];
XX = [xa ones(size(xa))];
% parameter estimation
theta = (X'*X+ eye(size(X,2))/sigma.^2)\(X'*Y);
% plotting
plot(xa,XX*theta);

%% MAXLIK using log basis
X  = [log(d_w) ones(size(d_h))];
XX = [log(xa) ones(size(xa))];
% parameter estimation
theta = (X'*X+ eye(size(X,2))/sigma.^2)\(X'*Y);
% plotting
plot(xa,XX*theta,'m--');
% writeout
datout('ml_log.res',[xa XX*theta],outFormat);

%% MAXLIK using fifth order pol basis
X  = [d_w.^5 d_w.^4 d_w.^3 d_w.^2 d_w ones(size(d_h))];
XX = [xa.^5  xa.^4  xa.^3  xa.^2  xa  ones(size(xa))];
% parameter estimation
theta = (X'*X+ eye(size(X,2))/sigma.^2)\(X'*Y);
% plotting
plot(xa,XX*theta,'k.-');
% writeout
datout('ml_five.res',[xa XX*theta],outFormat);

%% MAXLIK using ninth order pol basis
X  = [d_w.^9 d_w.^8 d_w.^7 d_w.^6 d_w.^5 d_w.^4 d_w.^3 d_w.^2 d_w ones(size(d_h))];
XX = [xa.^9  xa.^8  xa.^7  xa.^6  xa.^5  xa.^4  xa.^3  xa.^2  xa  ones(size(xa))];
% parameter estimation
theta = (X'*X+ eye(size(X,2))/sigma.^2)\(X'*Y);
% plotting
plot(xa,XX*theta,'c-.');
% writeout
datout('ml_nine.res',[xa XX*theta],outFormat,0,300);