library(gclus) # le package dans lequel se trouvent les données bank data(bank) # charge les données bank #---- 1a -------- cat('-----------------------------','\n') cat('On suppose les descripteurs de billets véritables & contrefaits distribués selon deux gaussiennes hétéroscédastiques','\n') cat('-----------------------------','\n') z<- bank[,1] x <- bank[,-1] n=nrow(x) # taille de l'échantillon n1=sum(z==0) # nbre de billets véritables n2=sum(z==1) # nbre de billets contrefaits pi1 <- n1/n # proportion a priori de billets véritables pi2 <- n2/n # proportion a priori de billets contrefaits xbar1 <- colMeans(x[z==0,]) # moyennes des descripteurs pour les billets véritables xbar2 <- colMeans(x[z==1,]) # moyennes des descripteurs pour les billets contrefaits W1 <- var(x[z==0,]) # variance des descripteurs pour les billets véritables W2 <- var(x[z==1,]) # variance des descripteurs pour les billets contrefaits Q <- function(x,pi,xbar,W){ # puisqu'on suppose les gaussiennes hétéroscédastiques, la fonction de score est quadratique as.matrix(x-xbar)%*%solve(W)%*%t(as.matrix(x-xbar))-2*log(pi)+log(det(W)) } #---- 1b -------- plot(x,col=c(rep('blue',n1),rep('red',n2)), main='classes réelles') cat('temps d arret pour observer les classes réelles ; appuyez sur une touche','\n') scan('') #---- 1c -------- zest_heteroscédastic <- matrix(nrow=n,ncol=1) scoretab <- matrix(0,nrow=n,ncol=2) for (i in 1:n){ scoretab[i,1]=Q(x[i,],pi1,xbar1,W1) scoretab[i,2]=Q(x[i,],pi2,xbar2,W2) } zest_heteroscédastic[scoretab[,1]=scoretab[,2]]=2 plot(x,col=zest_heteroscédastic, main='classes estimées (modèle gaussien hétéroscédastique)') cat('temps d arret pour observer les classes estimées (modèle gaussien hétéroscédastique) ; appuyez sur une touche','\n') scan('') #---- 1d -------- confusion <- table(z,zest_heteroscédastic) cat('tableau de confusion classe réelle x classe estimée','\n') print(confusion) erreur <- (confusion[1,2]+confusion[2,1])/n cat('erreur de classement (resubstitution) : ', erreur,'\n') # ---- 2. On suppose maintenant les descripteurs de billets véritables & contrefaits distribués selon deux gaussiennes homoscédastiques ---- #---- 2a -------- cat('-----------------------------','\n') cat('On suppose les descripteurs de billets véritables & contrefaits distribués selon deux gaussiennes homoscédastiques','\n') cat('-----------------------------','\n') # Les etimateurs des probabilités a priori ne changent pas pi1 <- n1/n # proportion a priori de billets véritables pi2 <- n2/n # proportion a priori de billets contrefaits # Les etimateurs des moyennes ne changent pas xbar1 <- colMeans(x[z==0,]) # moyennes des descripteurs pour les billets véritables xbar2 <- colMeans(x[z==1,]) # moyennes des descripteurs pour les billets contrefaits # Les etimateurs des variances changent W1 <- var(x[z==0,]) # variance des descripteurs pour les billets véritables W2 <- var(x[z==1,]) # variance des descripteurs pour les billets contrefaits W <- (n1*W1+n2*W2)/n W1 <- W2 <- W L <- function(x,pi,xbar,W){ # puisqu'on suppose les gaussiennes homoscédastiques, la fonction de score est linéaire as.matrix(x)%*%solve(W)%*%as.matrix(xbar)-xbar%*%solve(W)%*%as.matrix(xbar)/2+log(pi) } #---- 2b -------- plot(x,col=c(rep('blue',n1),rep('red',n2)), main='classes réelles') cat('temps d arret pour observer les classes réelles ; appuyez sur une touche','\n') scan('') #---- 2c -------- zest_homoscedastic <- matrix(nrow=n,ncol=1) scoretab <- matrix(0,nrow=n,ncol=2) for (i in 1:n){ scoretab[i,1]=L(x[i,],pi1,xbar1,W1) scoretab[i,2]=L(x[i,],pi2,xbar2,W2) } zest_homoscedastic[scoretab[,1]>=scoretab[,2]]=1 zest_homoscedastic[scoretab[,1]