%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% non_linear_fit.m
% An Octave sript to compute an non-linear
% fit to a set of data.
% Originally from work done in a numerical
% computing course this past summer.
% At the time, I couldn't solve the problem
% in FORTRAN, which I was using for the course.
%
% In general, the problem we want to solve is A dot a = b
% where the matrix A is defined as Aij=Xj(xi)/sigmai and 
% bi = yi/sigmai.  a is the coefficents of X which solve
% equation
%
% Daniel Garcia-Galili - dgg@marlboro.edu
% For Programming WOrkshop, Spring 2006
%
% THIS DOES NOT WORK RIGHT NOW!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



x = [1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 7.00, 8.00, 9.00, 10.00]*pi/5;

%Just cause we're going to have to count this high a lot
rows=length(x)+1;

%here's where we build our random data set
i=1;
while (i < rows)
  if (0.5 < rand)
    y(i) = sin(x(i))+.1*rand;
  else
    y(i) = sin(x(i))-.1*rand;
  endif
  i++;
endwhile

%And our sigma
sig = 0.1*rand(rows,1)

i=1;
while (i<rows)
  j=1;
  while (j<(rows)/2)
    A(i,j)=(x(i)^(j-1))/sig(i)
    j++;
  endwhile
  i++;
endwhile

i=1;
while (i<rows)
  b(i)=y(i)/sig(i);
  i++;
endwhile;

aa=(A'*A)*(A'*b);

%for i=(1:rows)
%  for j=(1:length(aa))
%    yfit(i)=aa(j)*x(i)^(j-1)
%  endfor
%endfor

i=1;
while (i<rows)
  yfit(i)=(aa(1)+(aa(2)*x(i))+(aa(3)*x(i)^2)+(aa(4)*x(i)^3)+(aa(5)*x(i)^4));
    i++;
endwhile

i=1;
while (i<rows)
  data(i,1)=x(i);
data(i,2)=y(i);
fit(i,1)=x(i);
fit(i,2)=yfit(i);
i++;
endwhile
data
fit
clearplot
gplot[0:2*pi][-1.25:1.25](data)
hold on
gplot(fit)
