Brique BASTA´Etude de cas :´Etude de la distribution des sauts d’amplitude dans une imagenum´eriqueAudran Le Baron17 d´ecembre 2002Fig. 1 – Image de travailLes donn´eesNoustravailleronstoutaulongdecette´etudedecassurl’imageaero.tiff,image en 256 niveaux de gris de taille 510×341. Nous chargeons cette imageet l’affichons par le code MatLab suivant. L’image est repr´esent´ee figure 1.a=imread(’aero.tiff’);imagesc(a);colormap gray;colorbar;axis equal;Nous travaillerons plus particuli`erement sur le vecteur des sauts d’ampli-tude horizontaux, not´e x, obtenu par le code suivant :a = double(a);[M, N] = size(a);h = a(:,2:end) - a(:,1:end-1); % Matrice des sauts horizontaux% (taille M x N-1)x = reshape(h, prod(size(h)), 1); % Matrice aplatie en vecteurLa figure 2 repr´esente les sauts d’amplitude horizontaux de cette image.1−50150 050 100 10050 1500 200−50 250300−100350−15040050 100 150 200 250 300 350 400 450 500Fig. 2 – Sauts d’amplitude horizontaux1 Mod´elisation par une gaussienneConsid´erons pour commencer le mod`ele param´etrique dans lequel la den-sit´edenotrevariableal´eatoire(valeurdusautd’amplitude)suituneloigaus-sienne : 21 (x−μ)g(x) = √ exp −22 2σ2πσ1.1 Comparaison d’histogrammesLe calcul de la moyenne et de la variance des observations se fait par lecode suivant :m = mean(x);v = var(x);Les r´esultats sont :μ = 0,01572σb = 337,3939Par le code suivant, on trace la distribution gaussienne dont les ...
Les donnees ´ Noustravailleronstoutaulongdecettee´tudedecassurl’image aero.tiff , image en 256 niveaux de gris de taille 510 × 341. Nous chargeons cette image etl’affichonsparlecodeMatLabsuivant.L’imageestrepre´sent´eefigure1. a=imread(’aero.tiff’); imagesc(a); colormap gray; colorbar; axis equal; Noustravailleronsplusparticulie`rementsurlevecteurdessautsd’ampli-tudehorizontaux,note´ x , obtenu par le code suivant : a = double(a); [M, N] = size(a); h = a(:,2:end) - a(:,1:end-1); % Matrice des sauts horizontaux % (taille M x N-1) x = reshape(h, prod(size(h)), 1); % Matrice aplatie en vecteur La fi re 2 r ´ te les sauts d’amplitude horizontaux de cette image. gu epresen
1
Fig. 2 – Sauts d’amplitude horizontaux
1Mod´elisationparunegaussienne Consid´eronspourcommencerlemod`eleparam´etriquedanslequelladen-site´denotrevariableale´atoire(valeurdusautd’amplitude)suituneloigaus-sienne : g ( x ) = √ 21 πσ 2 exp − ( x 2 − σ 2 µ ) 2 1.1 Comparaison d’histogrammes Le calcul de la moyenne et de la variance des observations se fait par le code suivant : m = mean(x); v = var(x); Le ´ ltats sont : s resu
µ = 0 , 0157 b σ 2 = 337 , 3939 Par le code suivant, on trace la distribution gaussienne dont les para-me`tressontlesmoyenneetvarianceempiriques,avecl’histogrammedesob-servations.Ler´esultatn’estpastr`essatisfaisantcommelemontrelafigure ?? .
2
Fig. 3 – Histogramme des observations et distribution gaussienne
x range = [-200:.1:200]; _ gauss_m_v = 1/(sqrt(2*pi*v)) * exp(-(x_range-m).^2 ./ (2*v)); hist(x, 100); hold on plot(x_range, length(x) * gauss_m_v, ’r’); hold off;
1.2 Comparaison de quantiles Soit F µ σ lafonctiondere´partitiond’unedistributiongaussiennedepa-, ram`etres( µ, σ ). Nous noterons F = F 1 , 1 . Le quantile d’ordre p de g µ,σ est le x tel que :
t − µ ) 2 p = F µ,σ ( x ) = Z − x ∞ √ 21 πσ 2 exp − (2 σ 2 dt 1 x − µ t 2 dt = Z −∞ σ exp − 2 σ = σ 1 F xσ − µ
3
Fig. 4 – Quantiles d’une loi gaussienne en fonction des quantiles d’une autre loi gaussienne
Ainsi, si q µ,σ ( p ) est le quantile d’ordre p d’une gausienne ( µ, σ ) et q ( p ) celuid’unegaussiennecentre´er´eduite,ona: q µ,σ ( p ) = σq ( σp ) + µ C’est pourquoi, lorsque l’on trace les quantiles d’une variable gaussienne par rapport aux quantiles d’une seconde variable gaussienne, on obtient une droite ( Cf. figure 4). Ce n’est pas le cas lorsque l’on trace les quantiles des observations ( Cf. figure 5) en fonction des quantile d’une loi normale par le code suivant, comme le montre la figure 6. g_m_v = m + sqrt(v) * randn(length(x), 1); qqnorm(g1); % C’est bien une droite
qqnorm(x); % Ce n’est PAS une droite On note qu’il y a trop de poids dans la queue de la gaussienne et pas assez au ✭ centre ✮ .
4
Fig. male
6
–
Quantiles
Fig. 5 – Quantiles des observations
des observations en fonction des quantiles
5
d’une
loi
nor-
2Mod´elisationparuneloi ✭ laplaciennege´-ne´ralisee ✮ ´ Changeonsdemode`leparam´etriquepourprendreceluid’uneloideLa-placege´n´eralise´e,dedensit´e: αη = h ( x )2Γ(1 /α )exp( − | ηx | α ) 2.1Estimationparlam´ethodedesmoments Unepremi`ereme´thoded’estimationestlam´ethodedesmoments.Nous avonsdeuxparam`etresscalaires( α et η ), donc deux moments qui sont ici la varianceetlekurtosis,donne´spar: σ 2 = η 2 ΓΓ 3 αα 1 Γ κ = α 1 Γ α 5 Γ 3 α 2 Le code suivant permet de tracer κ ( α ) et de calculer le kurtosis empirique κ . b alpha_vect = [25:100]/50; k = kurt(alpha_vect); plot(alpha_vect, k); k_emp = kurtosis(x); κ ( α )estrepre´sente´figure7etlecalculdukurtosisempiriquedonn b e κ = 17 , 8968.Doncd’apre`slacourbe, α se trouverait entre 0 , 5 et 0 , 6. D’o`ulare´solutionnume´riquesuivante: alpha_vect = [0.5:.0001:0.6]; [trash, indice] = min(abs(kurt(alpha_vect) - k_emp)) alpha = alpha_vect(indice) eta = 1 / sqrt(v * gamma(1/alpha) / gamma(3/alpha)) Lere´sultatdececalculdonne: α = 0 , 5672 η = 0 , 3469
Fig. 8 – Observations, gaussienne et Laplace (moments)
[temp, indice] = min(abs(repartition - p_range(i))); quant_emp(i) = x_range(indice); end plot(quantile(x, p_range), quant_emp); Onconstateunecertaineame´liorationparrapportaumode`legaussien (voir figure 9). Il y a cependant toujours trop de poids dans la queue de la distribution de Laplace.
2.2Estimationparlam´ethodedumaximumdevrai-semblance Uneautremethoded’estimationdesparam`etresestlam´ethodedumaxi-´ mumdevraisemblance.Ici,lesdeuxe´quationsdevraisemblancesont: n 1 /α + Ψ(1 /α ) /α 2 − n 1 X log( n | X i | ) | ηX i | α i =1 1 η = α i = n X 1 | X i | α ! − α n Parite´rationssuccessivesetgrˆaceaucodesuivant,oncalcule˜ α et η ˜.
8
Fig. 9 – Quantiles de la loi de Laplace en fonction des quantiles empiriques (moments)
% Suppression des zeros [trash1 trash2 y] = find(x);
% Initialisation n = length(x);
alpha = 0.5672; eta = 0.3469;
pas = 0.01; ecart = 0.5;
logvrais = -Inf;
% Boucle alpha_vect = [alpha - ecart : pas : alpha + ecart]; for i=1:30 [trash indice] = min(abs(eq4(y, alpha_vect, eta, n))); alpha = alpha_vect(indice); eta = power(alpha / n * sum (power(abs(x), alpha)), -1/alpha); newlogvrais = sum(log(laplace(x, alpha, eta)));
i alpha eta newlogvrais (newlogvrais - logvrais)/abs(newlogvrais) logvrais = newlogvrais; logvraisvect(i) = logvrais; end ; plot(logvraisvect); L’e´volutiondelalog-vraisemblaceaufuret`amesuredesit´erationsest ´nte´figure10.Aubutde18it´erations,lalog-vraisemblancesestabilise represe o pour les valeurs :
˜ 0 , 4072 α = ˜ 1 , 7715 η = Onpeutalors,pourcesnouvellesvaleursdesparam`etres,retracerladis-tribution et les quantiles. Ceci est fait aux figures 11, 12. On constat ` e a nouveauuneame´lioration,maislemode`len’esttoujourspasparfait.