%LEAST SQUARES ANALYSIS OF MAMMALIAN HEART RATE

%w body weights
%r corresponding heart rates

w = [3.5 4 6 25 103 117 200 252 300 437 1340 2000 2700 5000 22500 30000 33000 50000 70000 100000 415000 450000 500000 3000000];
r = [787 660 588 670 347 300 420 352 300 269 251 205 187 120 100 85 81 70 72 70 45 38 40 48];

x1 = w.^(-1/3)%raise each component of vector w to -1/3 power.
x2 = w.^(-2/3)

numerator = r*x1';%apply transpose operator ' to x1 to compute dot product.
denominator = sum(x2);%compute the sum of each component

k1 = numerator/denominator;

%%Now build the model on the first 2/3 of the data (16 points)
ws = w(8:24)
rs = r(8:24)

x1 = ws.^(-1/3)%raise each component of vector w to -1/3 power.
x2 = ws.^(-2/3)

numerator = rs*x1';%apply transpose operator ' to x1 to compute dot product.
denominator = sum(x2);%compute the sum of each component

k2 = numerator/denominator;%see formula in section 3.1.1

hold on
plot(w.^(-1/3),r,'o')%plot raw data
plot(w.^(-1/3),k1*w.^(-1/3),'--x')%plot first model
plot(w.^(-1/3),k2*w.^(-1/3),'--v')%plot second model

title('mammalian heart rate model')
xlabel('weight w^{(-1/3)}')
ylabel('pulse rate')
legend('raw data','least squares fit (all data)', 'least squares fit 2/3 data')