[heart_rate depth duration bird] ...
  = textread('Penguins.txt','%f%f%f%s','headerlines',1);

x = duration;
y = heart_rate;
b1 = sum((x-mean(x)).*(y-mean(y)))/sum((x-mean(x)).^2)
b1 = corr(x,y)*std(y)/std(x)
b0 = mean(y) - b1*mean(x)
	
figure(1);
plot(x, y, 'b.');
xlabel('Dive Duration (min)');
ylabel('Heart Rate (beats/min)');
title('Figure 3.19: Life as a Penguin');

hold on;
x0 = min(x);
x1 = max(x);
y0 = b0 + b1*x0;
y1 = b0 + b1*x1;
plot([x0 x1], [y0 y1], 'b-');
hold off;

yhat = b0 + b1*x;
r = y - yhat;
figure(2);
plot(yhat, r, 'b.');
xlabel('Predicted Values (y_hat)');
ylabel('Residuals');
title('Figure 3.20');

n = length(x);
A = [ ones(n,1), yhat, yhat.^2, yhat.^3, yhat.^4, yhat.^5 ];
b = A \ r;

hold on;
yhat = sort(yhat);
A = [ ones(n,1), yhat, yhat.^2, yhat.^3, yhat.^4, yhat.^5 ];
plot(yhat,A*b,'b-');
hold off;