newton.m

optimisation中的牛顿法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         TP1 Optimisation et commande                                     %
%          Recherche d'optimum--------Methode Newton + Moindres Carrees                    %
%          Programmation par   Nicolas MARTIN  et  Ning QU                                 %
%          Date: 11/10/2007                                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
%methode de Newton et moindre carree
%chargement des donnees
load('granul');

%valeur initiale
%estimation
a1=10;a2=8;
mu1=450;mu2=698;sig1=70;sig2=74;
teta=[mu1;mu2;sig1;sig2];
epsilon=0.0001;
k=2;
G=[2;2;2;2];
i=0;

% tableau d'evolution
evo_teta=[];
evo_T=[];

%boucle d'iteration
while norm(G)>epsilon
    i=i+1; %incremente le nombre d'iteration
    % Calcul du signal
    ym=a1*exp(-(x-mu1).^2/sig1^2)+ a2*exp(-(x-mu2).^2/sig2^2);
    e=y-ym; % erreur
    % Calcul du Gradient
    G1=-a1/sig1^2*e'*((x-mu1).*exp(-(x-mu1).^2/sig1^2));
    G2=-a2/sig2^2*e'*((x-mu2).*exp(-(x-mu2).^2/sig2^2));
    G3=-a1/sig1^3*e'*((x-mu1).^2.*exp(-(x-mu1).^2/sig1^2));
    G4=-a2/sig2^3*e'*((x-mu2).^2.*exp(-(x-mu2).^2/sig2^2));

    % Calcul de la derive seconde
    H1=4*a1^2/sig1^4*(((x-mu1).^2)'*exp(-(x-mu1).^2/sig1^2).^2);
    H2=4*a2^2/sig2^4*(((x-mu2).^2)'*exp(-(x-mu2).^2/sig2^2).^2);
    H3=4*a1^2/sig1^6*((x-mu1).^4)'*exp(-(x-mu1).^2/sig1^2).^2;
    H4=4*a2^2/sig2^6*((x-mu2).^4)'*exp(-(x-mu2).^2/sig2^2).^2;

    G=[G1;G2;G3;G4]; % vecteur
    H=[H1 0 0 0;0 H2 0 0;0 0 H3 0;0 0 0 H4]; %matrice Hessien

    %Calcul des nouvelle valeur
    teta=teta - 0.5*inv(H)*G;
    mu1=teta(1);
    mu2=teta(2);
    sig1=teta(3);
    sig2=teta(4);
    evo_teta(i,:)=teta';

    %calcul valeur moindre carree
    E1=exp(-(x-teta(1)).^2/teta(3)^2);
    E2=exp(-(x-teta(2)).^2/teta(4)^2);
    X=[E1 E2];
    T= inv(X'*X)*X'*y;
    a1=T(1);
    a2=T(2);
    evo_T(i,:)=T';
end
teta,i,T
%signal final
ym=a1*exp(-(x-mu1).^2/sig1^2)+ a2*exp(-(x-mu2).^2/sig2^2);
%signal bruite estime
ymb=a1*exp(-(x-mu1).^2/sig1^2)+ a2*exp(-(x-mu2).^2/sig2^2)+rand(size(x))/20;
%courbe bruite
plot(x,y,'*')
hold on
%courbe estime
plot(x,ym,'r')
% Courbe bruite estime
hold on
plot(x,ymb,'g')