solutions des exercices de programmation
Exo1:
l’observation set Y = + N:où est une constante et N v N (0; 1):
H0 : = 0
= 2 avec P [ = 2] = 1=2
H1 :
= 2 avec P [ = 2] = 1=2
1) Le rapport de vraissemblance s’écrit:
1 p1 (y 2)2 1 p1 (y+2)2
2 2 exp( 2 ) + 2 2 exp( 2 )
(y) = 2
p1
2
exp( y2 )
2
e
= [e2y + e 2y
]
2
2
= e cosh(2y)
le LRT est:
H1
(y) = e 2
cosh(2y) ?
H0
H1
cosh(2y) ? e2
H0
La fonction cosh(:) étant monotone et paire, la règle de décision peut être
mise sous forme:
H1
jyj ?
H0
1 2
(où = 2 arg cosh( e )). La règle de décision devient alors:
décider H0 : y
décider H1 : ( 1 y )U ( y +1)
La probabilité de fause alarme est donnée par:
Z Z +1
1 y2 1 y2
PF = p exp( )dy + p exp( )dy
1 2 2 2 2
Z +1
1 y2
= 2 p exp( )dy = 2Q( )
2 2
PF = 2Q( )
donc pour une proba de fausse alarme donnée, on peut tirer le seuil par:
1 PF
=Q ( )
2
1
La probabilité de détection est:
Z Z Z
+1 Z
+1
1 1 (y 2)2 1 1 (y + 2)2 1 1 (y 2)2 1
Pd = p exp( )dy + p exp( )dy + p exp( )dy +
2 2 2 2 2 2 2 2 2 2
1 1
Z
+1 Z
+1
1 (y 2)2 1 (y + 2)2
= p exp( )dy + p exp( )dy
2 2 2 2
Pd = Q( 2) + Q( + 2)
2) Les courbes ROCs sont données par:
1 PF 1 PF
P d = Q(Q ( ) 2) + Q(Q ( ) + 2)
2 2
N.B: sous matlab
p on peut calculer Q(x) en utilisant la fonction erf c(x):ona
Q(x) = 21 erf c(x= 2) et pour = Q(x) on calcul x = Q 1 ( ) par x = sqrt(2)
erf cinv(2 ):
3) simulation:
La stratégie à suivre est la suivante:
1- générer n échantillons sous H0 et n échantillons sous H1 :
2- considerer un seuil et determiner empiriquement PF et Pd de la façon
suivante:
nbre d0 echantillons en valeur absolue sous H0 >
PF =
n
nbre d0 echantillons en valeur absolue sous H1
Pd =
n
3-répeter l’étape 2 pour plusieurs valeur de :
programme matlab:
clear all
close all
%****** partie théorique******
pf=linspace(0,1,100);
x=sqrt(2)*erfcinv(pf);
pd=0.5*(erfc((x-2)/sqrt(2))+erfc((x+2)/sqrt(2)));
…gure(1)
plot(pf,pd,’-.’)
hold on
% ** partie simulation ******
clear all
n=1000; % nbre d’échantillons
H0=randn(1,n); % n échantillons de l’hypothèse H0: loi normale centrée et
réduite
%mu=2*randsrc(1,n); % genere +2 ou -2 equiprobablement
2
% mu=4*(randint(1,n)-0.5);% une autre manière de generer +2 ou -2 equiprob-
ablement
%% voici une autre manière de générer +2 ou -2 équiprobablement
%u = random(’unif’, 0, 1, 1,n);
mu=rand(1,n);
indice=…nd(mu>=0.5);
mu(1:n)=-2;
mu(indice)=2;
H1=mu+H0; % n échantillons de l’hypothèse H1
…gure(2)
subplot(2,1,1)
hist(H0,100);
subplot(2,1,2)
hist(H1,100);
t(1)=0; % seuil initiale
r=1;
while (t(r)<=max(abs(H0)))
ii=…nd(abs(H1)>=t(r));
jj=…nd(abs(H0)>=t(r));
pd(r)=length(ii)/n; % proba de détection
pf(r)=length(jj)/n; % proba de fausse alarme
t(r+1)=t(r)+0.001;
r=r+1;
end
…gure(1)
plot(pf,pd,’r’)
xlabel(’proba de fausse alarme’)
ylabel(’proba de détection’)
%** seuil pour pf=0.3****
alpha=0.3;
gammathe=sqrt(2)*erfcinv(alpha) % théorique
indice=…nd(pf==0.3);
%indice=…nd((0.299<=pf)&(pf<=0.3009));
gammaemp=t(indice) % empirique
legend(’théorique’,’empirique’)
Exo 2:
H1 : Y = 1+N
H0 : Y = 0+N
on a: N v N (0; 4):
les pdfs conditinnelles sont:
2
1 (y 0)
fY =H0 (y=H0 ) = p exp( 2
2 2
3
2
1 (y 1)
fY =H1 (y=H1 ) = p exp( 2
2 2
On a vu au cours que pour ce problème, la proba de fausse alarme et de
détection sont données par:
0
PF = Q( )
1
Pd = Q( )
où le seuil est donné par:
2
ln( ) + 0 1
= +
1 0 2
P0 (C10 C00 )
avec =
(1 P0 )(C01 C11 )
après substitution, en posant d = 1 0
:
ln( ) d
PF = Q( + )
d 2
ln( ) d
Pd = Q( )
d 2
4
les risques conditionnelles sont alors:
R0 ( ) = C00 + (C10 C00 )PF
R1 ( ) = C01 + (C11 C01 )Pd
et le risque de Bayes s’écrit:
RB ( ) = P0 R0 ( ) + (1 P0 )R1 ( )
avec les données du problème on trouve:
70P0
=
75(1 P0 )
d = 1:5
programme matlab:
clear
P0=[0:0.001:1];
eta=(70*P0)./(75*(1-P0));
sigma=2;
mu0=-1.5;
mu1=1.5;
d=(mu1-mu0)/sigma;
alpha=log(eta)/d;
pf=0.5*erfc((alpha+d/2)/sqrt(2));
pd=0.5*erfc((alpha-d/2)/sqrt(2));
R0=0.1+0.7*pf;
R1=1-0.75*pd;
Rb=P0.*R0+(1-P0).*R1;
plot(P0,Rb,’.-’);
indice=…nd(Rb==max(Rb));
P0(indice)
Rfb=P0.*R0(indice)+(1-P0).*R1(indice);
hold on
plot(P0,Rfb,’–’)
indice1=…nd(P0==0.500);
Rfb05=P0.*R0(indice1)+(1-P0).*R1(indice1);
plot(P0,Rfb05,’:r’)
indice2=…nd(P0==0.800);
Rfb08=P0.*R0(indice2)+(1-P0).*R1(indice2);
plot(P0,Rfb08,’-.’)
5
6