%% Nelineární metoda nejmenších čtverců
% f=f(x,\beta)=f(x,\beta_1,\dots,\beta_M)
% Data: y_j=f(x_j,\beta)+\varepsilon_j,\ j=0,\dots,N
% Určení odhadu parametrů \beta_i: minimalizace S(\beta)=\sum(y_j-f(x_j,\beta))^2
% Newtonův-Gaussův algoritmus:
% \beta^0 - počáteční odhad, \beta^{k+1}=\beta^k+\Delta\beta, \Delta y=y-f(x,\beta^k),
% Jacobiova matice J=(J_{j,i}),\ J_{j,i}=\frac{\partial f(x_j,\beta^k)}{\partial \beta_I}
% rovnice pro \Delta\beta: J^TJ\Delta\beta=J^t\Delta y
%% Příklad: Michalisův-Mentenové model kinematiky enzymů f(x,\beta)=\frac{\beta_1x}{\beta_2+x}
% Funkce pro model:

f=@(x,beta)beta(1)*x(:)./(beta(2)+x(:)); % vysledek je ve sloupci

% Simulovaná data pro \beta_1=3,\ \beta_2=1:

beta=[3;1];
N=20;
x=linspace(1,2,N+1);
fx=f(x,beta);
y=fx+0.1*randn(size(x(:)));
plot(x,y,'r*',x,fx,'b')

%
% Funkce pro Jacobiovu matici

JF=@(x,beta)[,];

%
% Počáteční iterace:

b=[;];

% Spočítáme 10 iterací

for k=1:10
    dy=y(:)-f(x,b);
    J=JF(x,b);
    db=J\dy % inv(J'*J)*J'*dy
    b=b+db
    pause(1)
end
plot(x,y,'r*',x,fx,'b--',x,f(x,b))

%
% Větší počet dat::

N=100;
x=linspace(1,2,N+1);
fx=f(x,beta);
y=fx+0.1*randn(size(x(:)));
plot(x,y,'r*',x,fx,'b')
b=[2;2];

%
% Spočítáme 10 iterací

for k=1:10
    dy=y(:)-f(x,b);
    J=JF(x,b);
    db=J\dy
    b=b+db
    pause(1)
end
plot(x,y,'r*',x,fx,'b--',x,f(x,b))

%
% Linearizovaný model: