%Assignment 2 - Math 411 - Spring 97 
%Objective is to illustrate difference between   
% a linear regression on a semilog plot and a fit to 
% the original non-linear model. This does both 
% a linear regression on the semi-log data as well as
% a brute force best fit for a bounding box centered
% at the linear best fit.
%This is for equations like 
%      y = a (b^x) for constants a,b.
% This is the same as 
%     N(t) = N(0) exp(k t) for a constant k.
% k>0 gives exponential growth 
% k<0 gives exponential decay.
%First read in data from a file
%   load b:semilog.dat
% that contains arrays x(i) and y(i)
% or else create these data using something like
%   for i+1:10
%       x(i)=i;
%       y(i)=3*2^i;
%   end;
% and then modify the y(i) data values as you like.
%To run the code, read in a command file that is
%this file of commands for MATLAB
%by first telling MATLAB to use the b: drive
%   !b:
%then telling it to run lesson3b
%   nonlinreg
%where nonlinreg.m is this file on your 
%floppy in the b: drive
who,pause
!c:
%Now read in the even number of partitions to 
% use for the brute force non-linear fitting
np=input('give even number of partitions for brute force fit:    ')
% First plot the raw data
plot(x,y,'+w')
title('x versus y'),pause
% Now plot the data on semilog scale
semilogy(x,y,'+')
title('x versus y semilog scale'),pause
% Now do a linear regression of the
% x versus log(y) plot
ly=log10(y)
c=polyfit(x,ly,1),pause
% here c(1) is the slope of the regression
% which for the above equation is log10(b)
% and c(2) is the intercept of the regression
% which for the above is log10(a)
% Now plot the regression line on plot
% of x versus log(y)
m=[min(x) max(x)];
ym=polyval(c,m);
plot(x,ly,'+',m,ym)
title('x versus log(y)'),pause
% Now to show original data and the fitted
% regression model, first set up the fitted
% model y values (called yfit)
len=length(x)
a=10^c(2)
b=10^c(1),pause
linr2=0.;
%Compute square error in linr2
for i=1:len
   yfit(i)=a*b^x(i);    
   linr2=linr2+(yfit(i)-y(i))^2;
end
linr2,pause
% Plot the data and the regression fit on
% original data scale
plot(x,y,'+w',x,yfit)
title('x versus y and fit'),pause
% Set up the number of partitions, here 
% np2 is half the number of partitions plus 1
np2=1+np/2;
% Put the nonlinear sum of the square errors
% in the variable nonlinr2
nonlinr2=10^12;
% Here tsa and tsb will contain the best current
% estimates of a and b from the nonlinear search
tsa=a;
tsb=b;
% In the box centered at the best fit a and b 
% from the linear regression, set up np partitions
% in a,b space so proceed to do (np+1)^2 calculations
% of the error and keep the lowest value in nonlinr2
for i=1:np+1
    ta=a+(i-np2)*.1*a;
for j=1:np+1
    tb=b+(j-np2)*.1*b;
      sum=0.;
      for k=1:len
          yfit1(k)=ta*tb^x(k);
          sum=sum+(yfit1(k)-y(k))^2;
          end;
      if sum < nonlinr2
          nonlinr2=sum;
          tsa=ta;
          tsb=tb;
        else;
        end;
end;
end;
% Print out the best fit values from the nonlinear
% regression and the sum of the square errors
tsa
tsb, pause
nonlinr2, pause
% To graph the best fit, put the fitted values
% in the array yfit1
for k=1:len
   yfit1(k)=tsa*tsb^x(k);
   end;
plot(x,y,'+w',x,yfit,x,yfit1)
title('x versus y and both fits'),pause
z(1)=log10(yfit1(1));
z(2)=log10(yfit1(len));
plot(x,ly,'+',m,ym,m,z)
title('x versus log(y) and fits'),pause