% Copyright (C) 1992-1994 Richard Shrager
% Copyright (C) 1992-1994 Arthur Jutan
% Copyright (C) 1992-1994 Ray Muzic
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% leasqrdemo
%
% An example showing how to use non-linear least squares to fit 
% simulated data to the function:
%
%      y = a e^{-bx}

% 2001-02-05 Paul Kienzle
%   * collected example into a single script

function leasqrdemo
  % generate test data
  t = [1:10:100]';
  p = [1; 0.1];
  data = leasqrfunc (t, p);

  rnd = [   0.352509
    -0.040607
    -1.867061
    -1.561283
     1.473191
     0.580767
     0.841805
     1.632203
    -0.179254
     0.345208 ];

  % add noise
  % wt1 = 1 /sqrt of variances of data
  % 1 / wt1 = sqrt of var = standard deviation
  wt1 = (1 + 0 * t) ./ sqrt (data); 
  data = data + 0.05 * rnd ./ wt1; 

  % Note by Thomas Walter <walter@pctc.chemie.uni-erlangen.de>:
  %
  % Using a step size of 1 to calculate the derivative is WRONG !!!!
  % See numerical mathbooks why.
  % A derivative calculated from central differences need: s 
  %     step = 0.001...1.0e-8
  % And onesided derivative needs:
  %     step = 1.0e-5...1.0e-8 and may be still wrong

  F = @leasqrfunc;
  dFdp = @leasqrdfdp; % exact derivative
  % dFdp = @dfdp;     % estimated derivative
  dp = [0.001; 0.001];
  pin = [.8; .05]; 
  stol=0.001; niter=50;
  minstep = [0.01; 0.01];
  maxstep = [0.8; 0.8];
  options = [minstep, maxstep];

  global verbose;
  verbose=1;
  [f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...
    leasqr (t, data, pin, F, stol, niter, wt1, dp, dFdp, options);

function y = leasqrfunc(x,p)
  % sprintf('called leasqrfunc(x,[%e %e]\n', p(1),p(2))
  % y = p(1)+p(2)*x;
  y=p(1)*exp(-p(2)*x);

function y = leasqrdfdp(x,f,p,dp,func)
  % y = [0*x+1, x];
  y= [exp(-p(2)*x), -p(1)*x.*exp(-p(2)*x)];