% a minimal test of least squares fitting in octave

SNR=100/3;  % 3% noise
a=1.2;      % "true" values
b=1.8;

x=linspace(0, 5); 
y = 1./(1 + a*x.^b) + randn(1,100)*0.05;

pkg load optim
% fit using leasqr (x, y, pin, F), where
%   x    - Vector (or matrix) of independent variables.
%   y    - Vector or matrix of observed values.
%   pin  - Vector of initial parameters to be adjusted by leasqr.
%   F    - Name of function or function handle. 
%          The function must be of the form y = f(x, p), with y, x, p of the form y, x, pin.
% Outputs: [f, p, cvg, iter, corp, covp, covr, stdresid, Z, r2] = leasqr (…)
%   f    - Column vector of values computed: f = F(x,p).
%   p    - Column vector trial or final parameters, i.e, the solution.
%   cvg  - Scalar: = 1 if convergence, = 0 otherwise.
%   iter - Scalar number of iterations used.
%   corp - Correlation matrix for parameters.
%   covp - Covariance matrix of the parameters.
%   ...
% A simple-minded interpretation of the results of the fit:
%   From http://originlab.com/www/helponline/Origin/en/UserGuide/The_Fit_Results.html
%   "The square root of a main diagonal value of the covariance matrix is the 
%    Standard Error of the corresponding parameter"

function y = f(x, p)
  y = 1./(1+p(1)*x.^p(2));
end

p_in = [0.5;1];
[yfit p_out cvg iter corp covp] = leasqr(x, y, p_in, "f");

if (cvg == 1)  
  disp(sprintf("convergence achieved in %d iterations\n   initial    final      +/-",iter))
  er=sqrt(diag(covp));
  disp([p_in p_out er])
end

convergence achieved in 5 iterations
   initial    final      +/-
   0.500000   1.169004   0.047758
   1.000000   1.898204   0.065892

%plot --format svg -w 600 -h 400

figure(1); clf(1);
  hold on;
  plot(x,y,...
     sprintf(";input data, true a=%.2f, b=%.2f;b-",a,b));
  plot(x,f(x,p_in),...
     sprintf(";initial guess: a=%.2f, b=%.2f;g-",p_in(1),p_in(2)));
  plot(x,f(x,p_out),...
     sprintf(";best fit: a=%.2f +/- %.2f, b=%.2f +/- %.2f;r-",p_out(1),er(1),p_out(2),er(2)));
  dy=y-f(x,p_out);
  plot(x,dy,...
     sprintf(";misfit, total misfit^2=%.4f;c-",sum(dy.^2)));  % could also use sumsq(dy)

  xlabel("x");
  ylabel("y = [ 1 + ax^b ]^{-1}");
  hold off;