%
% FIGURE3.M
%

% start from a clean sheet

clear all

% this random seed gives a dataset with a light elephant and a heavy hippo

randn('seed', 4567890);

% load data

load mammals.dat 
load dinosaurs.dat

% generate heteroscedastic data using conventional model

x_mammal           = mammals(:,1) + mammals(:,2);
log10y_mammal      = -1.0743 + 2.7331*log10(x_mammal) ...
                   + randn(size(x_mammal))*0.1088;
y_mammal           = 10.^log10y_mammal;

log10x = 0.5:0.01:3;
log10y = -1.0743 + 2.7331*log10x;
log10u = log10y + 2*0.1088;
log10l = log10y - 2*0.1088;

x = 10.^log10x;
y = 10.^log10y;
u = 10.^log10u;
l = 10.^log10l;

opts = optimset('Display', 'iter');
w    = fminunc(@loss, log2([2 2]), opts, x_mammal, y_mammal);
a    = 2.^w(1);
b    = 2.^w(2);
Y_mammal = a*x_mammal.^b;

figure(1);
clf;
set(axes, 'FontSize', 14)
plot(x_mammal, y_mammal, 'r+', ...
     x,        a*x.^b,   'b-.', ...
     x,        y,        'r-', ...
     x,        u,        'g--', ...
     x,        l,        'g--');
axis([0 1000 -2.5e+6 25e+6]);
xlabel('Humerus + femur (mm)');
ylabel('Body mass (g)');
text(500,15e+6,sprintf('%.4f X^{%.4f}',a,b), 'FontSize', 14);
line([300 500], [0.40e+7 0.23e+7], 'Color', 'black');
text(200,5e+6,sprintf('0.0843 X^{2.7331}',a,b), 'FontSize', 14);
line([600 800], [1.45e+7 0.55e+7], 'Color', 'black');
text(500,-0.5e+6, 'elephant', 'FontSize', 14);
line([650 850], [-0.0e+7 0.55e+7], 'Color', 'black');
legend('mammals', 'nlr model', '\mu', '\mu \pm 2\sigma', 'Location', 'NorthWest');
drawnow;
print -depsc2 fig3a.eps

sigma = std(Y_mammal - y_mammal);

figure(2)
clf;
set(axes, 'FontSize', 14)
plot(Y_mammal, (y_mammal - Y_mammal)/sigma, 'r+');
axis([-1e+6 8e+6 -6 +6]);
xlabel('Predicted body mass (g)');
ylabel('Standardised residual');
print -depsc2 fig3b.eps

% bye bye...

