vb=readtable('bertalanffy.csv');
plot(vb.edad, vb.longitud,"or");

syms Linf t0 r;
preds = Linf*(1-exp(-r*(vb.edad-t0)));
sumaCuadrados = sum((vb.longitud-preds).^2);

dL=diff(sumaCuadrados,Linf);
dt0=diff(sumaCuadrados,t0);
dr=diff(sumaCuadrados,r);

ndL = @(x) eval(subs(dL,{Linf,t0,r},{x(1),x(2),x(3)}))
ndt0 = @(x) eval(subs(dt0,{Linf,t0,r},{x(1),x(2),x(3)}))
ndr = @(x) eval(subs(dr,{Linf,t0,r},{x(1),x(2),x(3)}))
sistemaEc = @(x) [ndL(x), ndt0(x), ndr(x)]

% Resolvemos
sol=fsolve(sistemaEc,[340,-1.7,0.33])
%Comprobamos
sistemaEc(sol)
% Hacemos gráfica
plot(vb.edad,vb.longitud,"ro")
x=linspace(-3,14,200);
ajusteVB=@(x) sol(1)*(1-exp(-sol(3)*(x-sol(2))));   
fit=ajusteVB(x);
hold on
plot(x,fit,"-b")
yline(0);
xline(0);
hold off


% Arrancamos de otra aproximación inicial
sol=fsolve(sistemaEc,[340,-2.5,0.5])
sistemaEc(sol)
% Hacemos gráfica
plot(vb.edad,vb.longitud,"ro")
x=linspace(-3,14,200);
ajusteVB=@(x) sol(1)*(1-exp(-sol(3)*(x-sol(2))));   
fit=ajusteVB(x);
hold on
plot(x,fit,"-b")
yline(0);
xline(0);
hold off
% Solución inadecuada

% Cálculo de una aproximación inicial adecuada
y=vb.longitud(2:end)
x=vb.longitud(1:(end-1));
rectaRegresion=polyfit(x,y,1);
rectaRegresion
rapprox=-log(rectaRegresion(1))
plot(x,y,"or")
xr=linspace(200,350,200);
yr=rectaRegresion(1)*xr+rectaRegresion(2);
hold on
plot(xr,yr,"-b")
hold off;

% En realidad la recta anterior no se ha ajustado a las parejas
% de puntos [L(t),L(t+1)]. Para ello sería preciso interpolar:
vb=sortrows(vb);
x=vb.edad;
y=vb.longitud;
% Puede haber varios valores de y para un mismo x, en cuyo caso no se
% podría interpolar; detectamos los valores de x que se repiten y les 
% asignamos el promedio de todas las y que les corresponden:
xrep=unique(x(find(diff(x)==0)+1));
for k=1:length(xrep)
    xpos=find(x==xrep(k));
    yrep(k)=mean(y(xpos));
    x(xpos)=[];
    y(xpos)=[];
end
x=[x; xrep];
y=[y; yrep'];
% Ahora hacemos la interpolación:
ini=ceil(min(vb.edad));
fin=floor(max(vb.edad));
xq=ini:fin;
yq=interp1(x,y,xq);

x=yq(1:(end-1))
y=yq(2:end)
rectaRegresion=polyfit(x,y,1);
rectaRegresion
rapprox=-log(rectaRegresion(1))
plot(x,y,"or")
xr=linspace(200,350,200);
yr=rectaRegresion(1)*xr+rectaRegresion(2);
hold on
plot(xr,yr,"-b")
hold off;

plot(vb.edad,vb.longitud,"or",xq,yq,"xb")


% Calculamos una aproximación para t0 "a ojo"
t0aprox=5+(1/rapprox)*log(1-310/340)

% Resolvemos:
sol=fsolve(sistemaEc,[340,t0aprox,rapprox])
sistemaEc(sol)


% Solución usando lsqcurvefit (Ajuste de curvas por minimos cuadrados,
% "least squares" en inglés. Sigue siendo necesaria una aproximación
% inicial

fun=@(pars,xdata) pars(1)*(1-exp(-pars(3)*(xdata-pars(2))));
initial=[340,t0aprox,rapprox]
paramsAjustados = lsqcurvefit(fun,initial,vb.edad,vb.longitud)

t = linspace(min(vb.edad),max(vb.edad));
plot(vb.edad,vb.longitud,'ko',t,fun(paramsAjustados,t),'b-')
legend('Datos','Von Bertalanffy','Location','southeast')
title('Ajuste del modelo de Von Bertalanffy')
