% 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 tenth order pol basis
X  = [d_w.^10 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.^10 xa.^9  xa.^8  xa.^7  xa.^6  xa.^5  xa.^4  xa.^3  xa.^2  xa  ones(size(xa))];
for sigma=10.^[6:-1:0];
  % parameter estimation
  sigP = eye(size(X,2))/sigma.^2;
  sigP(end,end) = 10^10;
  theta = (X'*X+ sigP)\(X'*Y);
  % plotting
  plot(xa,XX*theta,'-');
  % writeout
  fname = ['map_ten_' num2str(round(log10(sigma))) '.res'];
  datout(fname,[xa XX*theta],outFormat,0,400);
end;