Commit 85993d50 authored by Ulmer Louis's avatar Ulmer Louis
Browse files

maj IML

parent 3fa881d0
No preview for this file type
No preview for this file type
function [xapp, yapp, xtest, ytest, indices] = SepareDataNfoldCV(x, y, Nfold, NumFold);
%% Generer les donnees app et test pour Nfold-CV
%% Valable pour N classes
%% x : donnees
%% y : labels
%% Nfold : nombre de decoupage
%% NumFold : numero de la portion a prendre pour les donnees de test
%% G2 - Fevrier 2006
if NumFold > Nfold
error('Num Portion sup a Nfold')
end
classcode = unique(y);
Nbclasse = length(classcode);
IndApp = [];
IndTest = [];
for i=1:Nbclasse
Indclass_i = find(y==classcode(i));
N_i = length(Indclass_i);
LongPortion_i = round(N_i/Nfold); %% longueur de chaque portion
IndTestClass_i = Indclass_i( LongPortion_i*(NumFold-1) + 1 : min(LongPortion_i*NumFold, N_i) );
IndTest = [IndTest; IndTestClass_i];
IndAppClass_i = setdiff(Indclass_i, IndTestClass_i);
IndApp = [IndApp; IndAppClass_i];
end
xapp = x(IndApp, :); yapp = y(IndApp);
xtest = x(IndTest, :); ytest = y(IndTest);
indices.app = IndApp;
indices.test = IndTest;
%% Anale2016
clear all
load('Archive-DonneesExam/datatrain.mat')
load('Archive-DonneesExam/datatest.mat')
%% Tester méthode ACP sur ces données
[valprop, U, moy] = mypca(Xa);
Ct = projpca(Xa , moy, U);
%% Diagnostique
%Pourcentage de variance expliquee : 100*(sum_{k=1}^d lambda_k)/(sum_{k=1}^D lambda_k)
vpc=100*cumsum(valprop)/sum(valprop);
figure(1)
bar(vpc,.5); hold on;
plot(100*valprop/valprop(1),'r-o');
title('Valeurs propres et leur importance en pourcentage d''information')
legend('Variance cumulee (en %)','100*valprop/valprop(1)')
xlabel('Nombre d''axes')
grid on
hold off
% Je choisis la methode 3
% Graphiquement, 95% de la variance est expliquee pour d=135
% Cette valeur de d est alors la valeur appropriee
% projpca fonction permettant de projeter une matrice
% de donnees X de taille NxD sur les d premiers vecteurs propres.
d=find(vpc>=95,1); % nombre de composantes principales qu'on veut garder
disp(['Les d=' num2str(d) ' premieres composantes principales representent ' num2str(vpc(d)) '% de la variance expliquee'])
fprintf('\n')
% Pour avoir 95% de variance explique je choisi graphiquement d=75
%% Affichage 2D
clf
figure(4)
plot(Ct(Ya==3,1),Ct(Ya==3,2),'bo')
hold on
plot(Ct(Ya==8,1),Ct(Ya==8,2),'rd')
%% Representer la reconstruction d’un des caracteres sur les 10 premieres composantes principales.
%d=75;
%On selectionne le nombre de composantes factorielles voulues
P=U(:,1:d);
C=Ct(:,1:d);
Xhat= reconstructpca(C,P,moy);
%% Preparation des variables
Yt(Yt==3)=-1;
Yt(Yt==8)=1;
Ya(Ya==3)=-1;
Ya(Ya==8)=1;
%% Validation Croisee Finale
%[xapp, yapp, xval, yval] = splitdata(xapp, yapp, 2/3);
K=3 %nombre de folds
meilleur_C=zeros(K,1);
vectC=logspace(-2,2,9) %echelle logarithmique
for k=1:1:K
%decoupage des donnees entree en K folds
[xapp_fold,yapp_fold,xval_fold,yval_fold]=SepareDataNfoldCV(Xhat,Ya,K,k);
meanx=mean(xapp_fold);
stdx=std(xapp_fold);
[xapp_fold,xval_fold,~,~] = normalizemeanstd(xapp_fold,xval_fold,meanx,stdx);
% evaluation d'un modèle
precision=zeros(length(vectC),1);
for i=1:1:length(vectC)
[wapp,b]=monsvmclass(xapp_fold,yapp_fold,vectC(i));
ypred=monsvmval(xval_fold,wapp,b);
%f(x) et y on le meme signe ?
ind_bon=find(ypred.*yval_fold>0);
n_bon=length(ind_bon); %nombre de valeurs bien classifiées
precision(i)=n_bon/length(yval_fold)
end
[val,pos] = max(precision);
meilleur_C(k)=(vectC(pos)); % On selectione les meilleurs C de chaques folds
end
C_final=mean(meilleur_C) % On fait la moyenne des k meilleurs C
%% Perfomance sur l'ensemble de test
meanx=mean(Xhat);
stdx=std(Xhat);
[Xhat,Xt,~,~] = normalizemeanstd(Xhat,Xt,meanx,stdx);
[wapp,b]=monsvmclass(Xhat,Ya,C_final);
ypred_test=monsvmval(Xt,wapp,b);
ypred_app=monsvmval(Xa,wapp,b);
%% Calcul de l'erreur de test
%f(x) et y on le meme signe ?
ind_bon=find(ypred_test.*Yt>0);
n_bon=length(ind_bon);
accuracy=n_bon/length(Yt)
erreur_test=1-accuracy
%% mise en forme des datas pour le calcul de la matrice de confusion matlab
ypred_test(ypred_test>0)=1;
ypred_test(ypred_test<0)=-1;
ypred_app(ypred_app>0)=1;
ypred_app(ypred_app<0)=-1;
%% Erreur de classification et Matrices de confusion
erreur_d_apprentissage=mean(ypred_app~=Ya)
erreur_de_test=mean(ypred_test~=Yt)
% Trace des matrices de confusion
Matrice_confusion_apprentissage=confusionmat(Ya,ypred_app)
Matrice_confusion_test=confusionmat(Yt,ypred_test)
\ No newline at end of file
%% Anale2016
clear all
load('Archive-DonneesExam/datatrain.mat')
load('Archive-DonneesExam/datatest.mat')
%% Tester méthode ACP sur ces données
[valprop, U, moy] = mypca(Xa);
Ct = projpca(Xa , moy, U);
%% Diagnostique
%Pourcentage de variance expliquee : 100*(sum_{k=1}^d lambda_k)/(sum_{k=1}^D lambda_k)
vpc=100*cumsum(valprop)/sum(valprop);
figure(1)
bar(vpc,.5); hold on;
plot(100*valprop/valprop(1),'r-o');
title('Valeurs propres et leur importance en pourcentage d''information')
legend('Variance cumulee (en %)','100*valprop/valprop(1)')
xlabel('Nombre d''axes')
grid on
hold off
% Je choisis la methode 3
% Graphiquement, 95% de la variance est expliquee pour d=135
% Cette valeur de d est alors la valeur appropriee
% projpca fonction permettant de projeter une matrice
% de donnees X de taille NxD sur les d premiers vecteurs propres.
d=find(vpc>=95,1); % nombre de composantes principales qu'on veut garder
disp(['Les d=' num2str(d) ' premieres composantes principales representent ' num2str(vpc(d)) '% de la variance expliquee'])
fprintf('\n')
% Pour avoir 95% de variance explique je choisi graphiquement d=75
%% Affichage 2D
clf
figure(4)
plot(Ct(Ya==3,1),Ct(Ya==3,2),'bo')
hold on
plot(Ct(Ya==8,1),Ct(Ya==8,2),'rd')
%% Representer la reconstruction d’un des caracteres sur les 10 premieres composantes principales.
d=75;
%On selectionne le nombre de composantes factorielles voulues
P=U(:,1:d);
C=Ct(:,1:d);
Xhat= reconstructpca(C,P,moy);
%% Preparation des variables
Yt(Yt==3)=-1;
Yt(Yt==8)=1;
Ya(Ya==3)=-1;
Ya(Ya==8)=1;
%% Validation Croisee bis
%[xapp, yapp, xval, yval] = splitdata(xapp, yapp, 2/3);
K=3 %nombre de folds
meilleur_C=zeros(K,1);
vectC=logspace(-2,2,9) %echelle logarithmique
for k=1:1:K
%decoupage des donnees entree en K folds
[xapp_fold,yapp_fold,xval_fold,yval_fold]=SepareDataNfoldCV(Xhat,Ya,K,k);
meanx=mean(xapp_fold);
stdx=std(xapp_fold);
[xapp_fold,xval_fold,~,~] = normalizemeanstd(xapp_fold,xval_fold,meanx,stdx);
% evaluation d'un modèle
precision=zeros(length(vectC),1);
for i=1:1:length(vectC)
[wapp,b]=monsvmclass(xapp_fold,yapp_fold,vectC(i));
ypred=monsvmval(xval_fold,wapp,b);
%f(x) et y on le meme signe ?
ind_bon=find(ypred.*yval_fold>0);
n_bon=length(ind_bon); %nombre de valeurs bien classifiées
precision(i)=n_bon/length(yval_fold)
end
[val,pos] = max(precision);
meilleur_C(k)=(vectC(pos)); % On selectione les meilleurs C de chaques folds
end
C_final=mean(meilleur_C) % On fait la moyenne des k meilleurs C
%% Perfomance sur l'ensemble de test
[wapp,b]=monsvmclass(Xhat,Ya,C_final);
ypred_test=monsvmval(Xt,wapp,b);
ypred_app=monsvmval(Xa,wapp,b);
%%
%f(x) et y on le meme signe ?
ind_bon=find(ypred_test.*Yt>0);
n_bon=length(ind_bon);
accuracy=n_bon/length(Yt)
error=1-accuracy
%% mise en forme des datas pour le calcul de la matrice de confusion matlab
ypred_test(ypred_test>0)=1;
ypred_test(ypred_test<0)=-1;
%%
Matrice_confusion_apprentissage(Ya,ypred_app)
Matrice_confusion_test=confusionmat(Yt,ypred_test)
\ No newline at end of file
function [CM] = confusionmatrice(Yt,ypred)
[ind1,val]=find(ypred>0);
[ind_moins1,val]=find(ypred<0);
ypred(ind1)=1;
ypred(ind_moins1)=-1;
TP=numel(find(ypred(ind1)==Yt(ind1)));
TN=numel(find(ypred(ind_moins1)==Yt(ind_moins1)));
FP=numel(find(ypred(ind1)~=Yt(ind1)));
FN=numel(find(ypred(ind_moins1)~=Yt(ind_moins1)));
CM=[TP FP; FN TN]
end
function [J,lam] = cout(H,x,y,C,ind,c,A,b,lambda)
% [J,yx] = cout(H,x,b,C,indd,c,A,b,lambda);
[n,m] = size(H);
X = zeros(n,1);
posok = find(ind > 0);
posA = find(ind==0); % liste des contriantes saturees
posB = find(ind==-1); % liste des contriantes saturees
% keyboard
X(posok) = x;
X(posB) = C(posB); %% 03/01/2002
J = 0.5 *X'*H*X - c'*X;% + lambda'*(A'*X-b);
% 0 normalement
% keyboard
% lam = y'*X;
lam = 0; %? ?liminer
\ No newline at end of file
%%
clc
close all
clear
%% Lecture des donnees
data=load('Archive-DonneesExam/datatrain.mat');
X=data.Xa;
Y=data.Ya; % labels
disp(['Donnees : datatrain.mat'])
disp(['Dimensions X : ' num2str(size(X))])
disp(['Dimensions Y : ' num2str(size(Y))])
fprintf('\n')
modalitesY = unique(Y);
Y(Y==modalitesY(1))=-1; % 3
Y(Y==modalitesY(2))=1; % 8
modalitesY = unique(Y);
disp(['Les labels Y ont pour modalites ' num2str(modalitesY(1)) ' (3) et ' num2str(modalitesY(2)) ' (8)'])
fprintf('\n')
effectif=sum(Y==modalitesY(1));
disp(['Effectif pour la modalite ' num2str(modalitesY(1)) ' : ' num2str(effectif)])
effectif=sum(Y==modalitesY(2));
disp(['Effectif pour la modalite ' num2str(modalitesY(2)) ' : ' num2str(effectif)])
% Meme s'il y a plus de donnees ayant le label -1, aucune des deux classes
% est reellement sous representee.
disp(['Effectif total = ' num2str(size(X,1))])
fprintf('\n')
% On considere que l'effectif total est petit donc on fait de la validation
% croisee pour :
% la methode des SVM lineaires pour determiner l'hyperparametre C optimal
%% 1. Ecrire un programme Matlab qui realise les operations suivantes :
%% - calculer la matrice de projection P de l'ACP
[valprop, P, moy_X] = mypca(X);
%% - projeter les donnees dans un espace de dimension q < 784 en utilisant P
% Cours ACP Slide 19
% Choix de la dimension d du sous-espace :
% Methode 1 : validation croisee
% Methode 2 : detection d'un coude sur le graphique des valeurs propres
% Methode 3 : on choisit d de sorte qu'un pourcentage fixe (par exemple 95% de la
% variance soit expliquee)
% Pourcentage de variance expliquee : 100*(sum_{k=1}^d lambda_k)/(sum_{k=1}^D lambda_k)
vpc=100*cumsum(valprop)/sum(valprop);
figure
bar(vpc,.5); hold on;
plot(100*valprop/valprop(1),'r-o');
title('Valeurs propres et leur importance en pourcentage d''information')
legend('Variance cumulee (en %)','100*valprop/valprop(1)')
xlabel('Nombre d''axes')
grid on
hold off
% Je choisis la methode 3
% Graphiquement, 95% de la variance est expliquee pour d=135
% Cette valeur de d est alors la valeur appropriee
% projpca fonction permettant de projeter une matrice
% de donnees X de taille NxD sur les d premiers vecteurs propres.
d=find(vpc>=95,1); % nombre de composantes principales qu'on veut garder
disp(['Les d=' num2str(d) ' premieres composantes principales representent ' num2str(vpc(d)) '% de la variance expliquee'])
fprintf('\n')
%% - choisir la valeur appropriee de q avec 2 ? q ? 100
% Normalement la bonne valeur pour q est 134. Comme on veut une valeur
% entre 2 et 100, q=100 est la valeur la plus appropriee
d=100;
Ct = projpca(X, moy_X, P(:,1:d));
%% - apprendre un modèle SVM lineaire f (z) = w > z + b avec z la projection
%% du point x par la matrice P
C = logspace(-2, 2, 9); % choisi par moi, arbitraire
errMin = Inf;
bestC = [];
best_numEnsemble = [];
N = length(Y);
% On prepare les donnees pour la validation croisee :
% Cours Introduction au Machine Learning : Slides 27-28 : il est dit qu'on
% separe generalement les donnees en 5 pour la validation croisee.
% Cours Selection et evaluation de modeles : Slide 23
% A : Apprentissage
% V : Validation
% 1 ... N/5 ... 2N/5 ... 3N/5 ... 4N/5 ... N
% V A A A A
% A V A A A
% A A V A A
% A A A V A
% A A A A V
% ATTENTION ce code est un peu different du code des TP10 et TP12
nbEnsembles=5;
for numEnsemble = 0:nbEnsembles-1
intervalle = 1+numEnsemble*floor(N/nbEnsembles):(numEnsemble+1)*floor(N/nbEnsembles);
coord = logical(zeros(length(Y),1));
coord(intervalle)=1;
Cval = Ct(coord,:);
yval = Y(coord);
Capp = Ct(~coord,:);
yapp = Y(~coord);
% Normalisation des donnees
m_app=mean(Capp);
sigma_app=std(Capp);
[Capp_normalise,Cval_normalise,~,~] = normalizemeanstd(Capp,Cval,m_app,sigma_app);
% [xapp,xtest,meanxapp,stdxapp] = normalizemeanstd(xapp,xtest,meanx,stdx)
disp(['Donnees xval (numEnsemble=' num2str(numEnsemble) ') :'])
% En CM, en faisant une annale, on a dit :
% Afin de choisir le parametre C du SVM, il faut evaluer chaque modele
% (avec des C differents) sur les donnees de validation
% (pas sur les donnees d'apprentissage, pas sur les donnees de test)
for i = 1:9
[omega,b,alpha] = monsvmclass(Capp_normalise,yapp,C(i));
% [w, b, alpha] = monsvmclass(X, labels, C)
ypredit_val = monsvmval(Cval_normalise, omega, b);
% ypred = monsvmval(Xv, w, b)
% Erreur de validation
erreur_val(i,numEnsemble+1)=mean(ypredit_val~=yval);
disp(['erreur_val=' num2str(erreur_val(i,numEnsemble+1)) ' ; C=' num2str(C(i))])
if (erreur_val(i,numEnsemble+1) < errMin)
errMin = erreur_val(i,numEnsemble+1);
bestC = C(i);
best_numEnsemble = numEnsemble;
end
end
fprintf('\n')
end
disp(['Le meilleur C est C=' num2str(bestC) ' ; numEnsemble=' num2str(best_numEnsemble) ' ; errMin_val=' num2str(errMin)])
figure
hold on
grid on
for numEnsemble = 0:nbEnsembles-1
plot(C', erreur_val(:,numEnsemble+1), '-o')
end
title('Erreur de validation en fonction de C')
legend('numEnsemble=0','numEnsemble=1','numEnsemble=2','numEnsemble=3','numEnsemble=4')
hold off
%% 2. Evaluer le modèle sur les donnees d'apprentissage et de test en
%% affichant les matrices de confusion et les erreurs de classification
%% respectives. Les donnees de test sont disponibles dans le fichier
%% datatest.mat sous la forme X t ? R N t ×784 , Y t ? R Nt et chaque
%% etiquette yj ? {3, 8}.
intervalle = 1+best_numEnsemble*floor(N/nbEnsembles):(best_numEnsemble+1)*floor(N/nbEnsembles);
coord = logical(zeros(length(Y),1)); % ATTENTION length(Y) n'est plus egale a N
coord(intervalle)=1;
Cval = Ct(coord,:);
yval = Y(coord);
Capp = Ct(~coord,:);
yapp = Y(~coord);
% Normalisation des donnees
m_app=mean(Capp);
sigma_app=std(Capp);
[Capp_normalise,Cval_normalise,~,~] = normalizemeanstd(Capp,Cval,m_app,sigma_app);
% [xapp,xtest,meanxapp,stdxapp] = normalizemeanstd(xapp,xtest,meanx,stdx)
[omega,b,alpha] = monsvmclass(Capp_normalise,yapp,bestC);
% [w, b, alpha] = monsvmclass(X, labels, C)
ypredit_app = monsvmval(Capp_normalise, omega, b);
% ypred = monsvmval(Xv, w, b)
%% Matrice confusion : donnees d'apprentissage
%ypredit_app = monsvmval(Capp_normalise, omega, b);
% ypred = monsvmval(Xv, w, b)
classe1=1; % positifs
classe2=-1; % negatifs
[confusionMatrix_app, myCell_app, fpr_app, tpr_app, erreur_app]=matriceConfusion(yapp,ypredit_app,classe1,classe2);
%[confusionMatrix, myCell, fpr, tpr, erreur]=matriceConfusion(ytrue,ypred,classe1,classe2)
myCell_app
erreur_app
%%
C=projpca(X, moy_X, P(:,1:d));;
data = load('datatest.mat');
Xtest = data.Xt;
Ytest = data.Yt;
[valprop, P, ~] = mypca(Xtest);
moy_Xtest=mean(Xtest);
Ctest = projpca(Xtest, moy_Xtest, P(:,1:d));
modalitesY = unique(unique(Ytest));
Ytest(Ytest==modalitesY(1))=-1; % 3
Ytest(Ytest==modalitesY(2))=1; % 8
moy_C=mean(C);
sigma_C=std(C);
% Normalisation des donnees de tests
[~,Ctest_normalise,~,~] = normalizemeanstd(C,Ctest,moy_C,sigma_C);
ypredit_test = monsvmval(Ctest_normalise, omega, b);
erreur_test=mean(ypredit_test~=Ytest);
[confusionMatrix_test, myCell_test, fpr_test, tpr_test, erreur_test]=matriceConfusion(Ytest,ypredit_test,classe1,classe2);
% [confusionMatrix, myCell, fpr, tpr, erreur]=matriceConfusion(ytrue,ypred,classe1,classe2)
myCell_test
erreur_test
\ No newline at end of file
function [confusionMatrix, myCell, fpr, tpr, erreur]=matriceConfusion(ytrue,ypred,classe1,classe2)
modalites=[classe1 classe2];
n=length(modalites);
confusionMatrix=zeros(n,n);
for i=1:n
for j=1:n
confusionMatrix(i,j)=sum((ytrue==modalites(j)).*(ypred == modalites(i)));
end
end
TP=confusionMatrix(1,1);
FN=confusionMatrix(2,1);
FP=confusionMatrix(1,2);
TN=confusionMatrix(2,2);
fpr=FP/(FP+TN);
tpr=TP/(TP+FN);
% Formule trouvee sur :
% http://pageperso.lif.univ-mrs.fr/~remi.eyraud/CAD/theorie.pdf
% Ctrl Matrice de confusion
erreur = (FP+FN)/(FP+FN+TP+TN);
myCell=cell(4,3);
myCell{1,1} = 'Predite y^ \ Reelle y';
myCell{2,1} = 'C1 = Pos';
myCell{3,1} = 'C2 = Neg';
myCell{4,1} = 'Total';
myCell{1,2} = 'C1 = Pos';
myCell{1,3} = 'C2 = Neg';
for i = 1:n
for j = 1:n
myCell{i+1,j+1} = confusionMatrix(i,j);
end
end
myCell{4,2} = sum(confusionMatrix(:,1));
myCell{4,3} = sum(confusionMatrix(:,2));
end
\ No newline at end of file
function [xnew, lambda, pos, mu] = monqp(H,c,A,b,C,l,verbose,xinit)
% function [xnew, lambda, pos, mu] = monqp(H,c,A,b,C,l,verbose,X,ps,xinit)
%
% min 1/2 x' H x - c' x
% x
% contrainte A' x = b
%
% et 0 <= x <= C
%
% ENTREES
% x : vecteur de taille n qu'on cherche
% H : matrice de dimensions n x n
% c : vecteur de taille n
% A : matrice de taille m x n avec m : nombre de contraintes ?galit?s
% b : vecteur de taille m
% C : vecteur de taille n tq 0 <= x_i <= C_i avec C_i le ieme ?l?ment de C
% l : param?tre de stabilisation (joue le m?me role que lambda dans la
% : r?gression logistique)
% verbose : = 1 - on affiche les it?rations ; 0 - on affiche rien
% xinit : initialisation de x. Par d?faut xinit = zeros(n,1);
%
% SORTIES
% xnew : solution qui ne contient que les ?l?ments xi non nuls
% lambda : param?tres de lagrange associ?s ? la contrainte ?galit? A'*x = b
% : pour un svm de type w'*x + b on a b = lambda
% pos : indique les index des ?l?ments non nuls de x
% : ?ad xnew = x(pos);
% mu : param?tres de lagrange associ?s aux contraintes 0 <= x <= C
%
% Cr?dits : J.P. Yvon (INSA de Rennes) pour l'algorithme
% O. Bodard (Clermont Ferrand) pour le debugage on line
%
% S CANU - scanu@insa-rouen.fr
% Mai 2001
%--------------------------------------------------------------------------
% verifications
%--------------------------------------------------------------------------
[n,d] = size(H);
[nl,nc] = size(c);
[nlc,ncc] = size(A);
[nlb,ncb] = size(b);
if d ~= n
error('H must be a squre matrix n by n');
end
if nl ~= n
error('H and c must have the same number of row');
end
if nlc ~= n
error('H and A must have the same number of row');
end
if nc ~= 1