% OBJECTIVE FUNCTION
f = @(x,b) b(1)./(1+exp(b(2)-b(3).*x)).^(1./(b(4)))

% PROBLEM SETUP

dim = 4 % dimension of the minimization problem
nb = 15 % number of points

% NIST certified values
target = zeros(dim,1);
target(1) = 6.9964151270E+02;
target(2) = 5.2771253025E+00;
target(3) = 7.5962938329E-01;
target(4) = 1.2792483859E+00;

% start 1
ig1 = [100;10;1;1]

% start 2
ig2 = [75;5;0.75;1.3]

% points to fit
x = zeros(nb,1)
y = zeros(nb,1)
x(1)=1;	y(1)=16.08
x(2)=2;	y(2)=33.83
x(3)=3;	y(3)=65.8
x(4)=4;	y(4)=97.2
x(5)=5;	y(5)=191.55
x(6)=6;	y(6)=326.2
x(7)=7;	y(7)=386.87
x(8)=8;	y(8)=520.53
x(9)=9;	y(9)=590.03
x(10)=10;	y(10)=651.92
x(11)=11;	y(11)=724.93
x(12)=12;	y(12)=699.56
x(13)=13;	y(13)=689.96
x(14)=14;	y(14)=637.56
x(15)=15;	y(15)=717.41


% GRAPHIC REPRESENTATION
scatter(x,y)

% OPTIMIZATION 
rng(0)	% for tractability 
[b,info] = lsqcurvefit(f,ig1,x,y)	% running minimization 

% ANIMATION
% Uncomment the following code to generate animation % 
% rng(0)
% [~,info] = lsqcurvefit(f,ig1,x,y)
% info.plot = 'fit'
% info.animate = true	% plot animation 
% optimview('lsqcurvefit',info) 