// lecture des données
donnees = read("donnees.txt", -1, 2)
Xexp = donnees(:, 1);
Yexp = donnees(:, 2);
// regression
[S, Aopt] = leastsq(list(res, Xexp, Yexp),Ainit)
// loi modèle
Ymod = f(Aopt(1), Aopt(2), Xexp);
// affichage
plot(Xexp, Ymod, "-r")
plot(Xexp, Yexp, "+b")
xstring(0, ordonnee, "a = "+string(Aopt(1))+" pour "+string(pente)...
+" ; b = "+string(Aopt(2))+" pour "+string(ordonnee)...
+" ; variance = "+string(S))
 Exemple 2
Nous travaillons ici avec un pic généré par la somme d'une fonction gaussienne et d'une fonction lorentzienne. Ce genre de profil est classiquement utilisé pour décrire la forme des pics obtenus par diffractométrie X. On fait ce que l'on appelle une désommation de pics (on parle souvent de « déconvolution », mais ce terme est abusif puisqu'il ne s'agit pas de produit de convolution).
clf()
// Constantes
paramlorentz(1) = 5; // hauteur de la courbe lorentzienne
paramlorentz(2) = 2; // largeur de la courbe lorentzienne
paramgauss(1) = 10; // hauteur de la courbe gaussienne
paramgauss(2) = 3; // largeur de la courbe gaussienne
var=0.5; // variance de la loi normale
nbpts = 100 // nombre de points
demielargeur = max(3*paramgauss(2), 3*paramlorentz(2)) // pour intervalle x
pas = 2*demielargeur/nbpts;
Ainit = [1;1;1;1]; // paramètres initiaux du modèle
// Initialisation
X = zeros(nbpts,1);
Y = zeros(nbpts,1);
Yopt = zeros(nbpts,1);
// Fonctions
// Fonction gaussienne centrée sur 0
function [y] = gauss(A, x)
 // A(1) : hauteur de pic
 // A(2) : "largeur" du pic
 y = A(1)*exp(-x.*x/A(2));
endfunction
// Fonction lorentzienne centrée sur 0
function [y] = lorentz(A, x)
 // A(1) : hauteur de pic
 // A(2) : "largeur" du pic
 foo = A(2)*A(2)
 y = foo*A(1)./(foo + 4*x.*x);
endfunction
function [y] = profil(A, x)
 // A(1) : hauteur de pic lorentzien
 // A(2) : "largeur" du pic lorentzien
 // A(3) : hauteur de pic gaussien
 // A(4) : "largeur" du pic gaussien
 L = A(1:2);
 G = A(3:4);
 y = lorentz(L, x) + gauss(G, x);
endfunction
// Résidus
function [e] = residus(A, x, y)
 e = profil(A, x) - y;
endfunction
// Nuage de points
i=(1:nbpts)';
X = pas*i - demielargeur;
Y = profil([paramlorentz;paramgauss], x) + var*rand(X, "normal");
// Programme principal
[S, Aopt] = leastsq(list(residus, X, Y), Ainit)
Yopt = profil(Aopt, X);
YLopt = lorentz(Aopt(1:2),X);
YGopt = gauss(Aopt(3:4),X);
// affichage
plot(X, Yopt, "-r")
plot(X, YLopt, "-c")
plot(X, YGopt, "-g")
plot(X, Y, "+b")
hauteur = paramlorentz(1)+paramgauss(1);
xstring(-demielargeur, hauteur...
 "lorentzienne : Al = "+string(0.01*round(100*Aopt(1)))...
+" ; Hl = "+string(0.01*round(100*Aopt(2))))
xstring(-demielargeur, 3*hauteur/4...
 "gaussienne : Ag = "+string(0.01*round(100*Aopt(3)))...
+" ; Hg = "+string(0.01*round(100*Aopt(4))))
xstring(-demielargeur, hauteur/2...
 "variance : S = "+string(0.01*round(100*S))))
L'algorithme par défaut de codice_16 est l'algorithme de pseudo-Newton. On peut aussi utiliser l'algorithme du gradient conjugué en ajoutant codice_22 aux paramètres, ou bien algorithme foncitonnant avec les fonctions non-différentiables avec codice_23. On peut enfin utiliser la commande codice_24, qui utilise l'algorithme de Levenberg-Marquardt.
Notes.
"Interpolation, extrapolation et lissage" < ↑ > "Calcul différentiel et intégral"

Goélette Cardabela/Références/Référencement automatique

Goélette Cardabela/Référencement semi-automatique