function [r, J]=randJac(x)
    f = exp(A*x);
    J = [f, t.*f, t.^2.*f];
    r = f-y;
end
http_get("https://moodle.utc.fr/mod/resource/view.php?id=109806","dataGaussian.sod",follow=%t);
load dataGaussian.sod

clf
plot(t,y,"o")

// régression utilisant le "log trick"

A = [ones(t) t t.^2];
x = A\log(y);

tt = linspace(-3,9,1000);
plot(tt,exp(x(1)+x(2)*tt+x(3)*tt.^2),"color",[0 0.7 0],"thickness",2)

// les valeurs obtenues ci-dessus sont utilisées comme valeurs
// initiales pour l'algorithme de Gauss-Newton.
// Ces valeurs sont suffisamment proches des valeurs optimales
// pour ne pas avoir à utiliser la méthode de Levenberg-Marquardt

ITMAX = 1000;
TOL = 1e-15;
for k=1:ITMAX
    [r,J] = randJac(x);
    dx = J \ r;
    x = x-dx;
    if norm(dx) < TOL
        break
    end
end
disp(k)

plot(tt,exp(x(1)+x(2)*tt+x(3)*tt.^2),"r","color",[1 0 0],"thickness",2)
legend("data","log trick","native")

// les valeurs des paramètres "naturels" ne sont pas utilisées 

a = exp(x(1)-x(2)^2/4/x(3));
mu = -x(2)/2/x(3);
sigma = sqrt(-1/x(3));