% LMex110220_4th.m
% using as a model
% LMcs110204_4th.m
% Larsen & Marx Case Study 11.2.4, an example of exponential regression
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% Written by Eugene.Gallagher@umb.edu; written 11/23/10; revised 11/23/10
% DATA=[0 4.5e3;3 2.9e4;7 9e4;10 2.29e5;14 1.2e6;18 3.1e6;20 5.5e6];
DATA=[1 1 2 2 2 3 5 6 6 7
    94.5 86.4 71.0 80.5 81.4 67.4 49.3 46.8 42.3 36.6]';
[r,c]=size(DATA);
X=[ones(r,1) DATA(:,1)];Y=log(DATA(:,2));
B=X\Y;
fprintf('\nThe Y intercept = %8.3f and slope = %5.3f\n',B(1),B(2));
s=sprintf('y=%6.1f e^{%5.3f x}',exp(B(1)),B(2));
Yest=X*B;
Resid=Y-Yest;
% sort by the independent variable so the lines are smooth
[Xs,k]=sort(X(:,2));
% sort the estimated values using the index vector k:
Ys=Yest(k);
plot(X(:,2),Y,'o',Xs,Ys,'m',...
                        'LineWidth',2,...
                       'MarkerEdgeColor','k',...
                       'MarkerFaceColor','b',...
                       'MarkerSize',8) 
xlabel('Days after Injection','Fontsize',20)
ylabel('Serum Gold % Concentration','Fontsize',20)
axis([0.7 7.3 3.55 4.6])
figure(gcf)
pause
% Generated a fine-scale estimated curve
xintrp=[0:0.1:20]';
[r,c]=size(xintrp);
xintrp=[ones(r,1) xintrp];
yestintrp=xintrp*B;
plot(DATA(:,1),DATA(:,2),'o',xintrp(:,2),exp(yestintrp),'--r',...
                        'LineWidth',2,...
                       'MarkerEdgeColor','k',...
                       'MarkerFaceColor','b',...
                       'MarkerSize',8) 
xlabel('Days after Injection','Fontsize',20),
ylabel('Serum Gold % Concentration','Fontsize',20)
text(5,3.5e6,s,'Fontsize',20);
title('Exercise 11.2.20','Fontsize',22)
axis([0 8 0 120])
figure(gcf)
pause
plot(Yest,Resid,'ok');
xlabel('Expected value','Fontsize',20)
ylabel('Residual','Fontsize',20);
grid;figure(gcf)

% Use the Matlab statistical toolbox's regress for confidence limits
[B,BINT] = regress(Y,X);
fprintf(...
'The half life for removal of gold is %4.2f days\n',...
abs(log(2)./B(2))); 
% What is the confidence limit?
fprintf(...
'The 95%% confidence limits for the half life is: %3.2f and %3.2f days\n',...
abs(fliplr(log(2)./BINT(2,:)))); 

% This can also be tested for lack of fit using Gallagher's leastsqu.m
[b,Yest,resid,SS,F,R2,Vb,VYest]=leastsqu(X(:,2),Y,0,0);
% Lack of fit was non-signficant.