# setwd('/media/alourme/GERMAIN/DESPEG/2015-16/S1/master/scoring/Rcode/Rcode1/') library(Rmixmod) # la librairie Rmixmod data(finance) # les données finance contenues dans Rmixmod names(finance) # pour connaître le nom des variables disponibles : deux variables qualitatives (année & santé) & quatre ratios financiers (variables quantitatives) traindata <- finance[finance$Year=='2002',-c(1,2)] # données d'apprentissage : entreprises de 2002 décrites par les quatre ratios financiers trainlabels <- finance[finance$Year=='2002',2] # la santé des entreprises de 2002 (= le label). testdata <- finance[finance$Year=='2003',-c(1,2)] # données de test : entreprises de 2003 décrites par les quatre ratios financiers testlabels <- finance[finance$Year=='2003',2] # la santé des entreprises de 2003 (la variable qu'il faudrait retrouver grâce au classifieur) # 'Gaussian_pk_L_C' désigne le modèle gaussien homoscédastique # 'Gaussian_pk_Lk_Ck' désigne le modèle gaussien hétéroscédastique # pour connaître les autres modèles disponibles, consulter les deux articles [12] et [13] cités en références bibliographiques # statistiques des modèles homoscédastiques et hétéroscédastiques learn_homo <- mixmodLearn(traindata, knownLabels=trainlabels, models = mixmodGaussianModel(listModels=c("Gaussian_pk_L_C")), # homoscédastique criterion = c("BIC","CV")) learn_hetero <- mixmodLearn(traindata, knownLabels=trainlabels, models = mixmodGaussianModel(listModels=c("Gaussian_pk_Lk_Ck")), # hétéroscédastique criterion = c("BIC","CV")) # comparaison selon BIC cat('---- comparaison selon BIC -----','\n') cat('BIC du modèle homoscedastique : ',learn_homo[8][3][1],'\n') cat('BIC du modèle hétéroscedastique : ',learn_hetero[8][3][1],'\n') # comparaison selon l'erreur estimée par resubstitution cat('---- comparaison selon l erreur de classement -----' ,'\n') prediction_homo <- mixmodPredict ( data = testdata , classificationRule = learn_homo["bestResult"]) # classification des entreprises de 2002 selon le modèle homoscédastique pred_homo <- as.matrix(prediction_homo[5]) confusion <- table(pred_homo,testlabels) error_rate_homo_resubstitution=(confusion[1,2]+confusion[2,1])/sum(confusion) cat('erreur de classement du modèle homoscédastique (resubstitution) : ',error_rate_homo_resubstitution,'\n') prediction_hetero <- mixmodPredict ( data = testdata , classificationRule = learn_hetero["bestResult"]) # classification des entreprises de 2002 selon le modèle hétéroscédastique pred_hetero <- as.matrix(prediction_hetero[5]) confusion <- table(pred_hetero,testlabels) error_rate_hetero_resubstitution=(confusion[1,2]+confusion[2,1])/sum(confusion) cat('erreur de classement du modèle hétéroscédastique (resubstitution) : ',error_rate_hetero_resubstitution,'\n') # comparaison selon l'erreur estimée par VC cat('VC du modèle homoscedastique : ',learn_homo[8][3][2],'\n') cat('VC du modèle hétéroscedastique : ',learn_hetero[8][3][2],'\n') # comparaison selon l'erreur de classement des entreprises de 2003 cat('---- comparaison selon l erreur de classement des entreprises de 2003 -----','\n') # classification des entreprises de 2003 selon le modèle homoscédastique prediction_homo <- mixmodPredict ( data = testdata , classificationRule = learn_homo["bestResult"]) pred_homo <- as.matrix(prediction_homo[5]) confusion <- table(pred_homo,testlabels) error_rate_homo=(confusion[1,2]+confusion[2,1])/sum(confusion) cat('erreur de classement des entreprises de 2003 par le modèle homoscédastique : ',error_rate_homo,'\n') # classification des entreprises de 2003 selon le modèle hétéroscédastique prediction_hetero <- mixmodPredict ( data = testdata , classificationRule = learn_hetero["bestResult"]) pred_hetero <- as.matrix(prediction_hetero[5]) confusion <- table(pred_hetero,testlabels) error_rate_hetero=(confusion[1,2]+confusion[2,1])/sum(confusion) cat('erreur de classement des entreprises de 2003 par le modèle hétéroscédastique : ',error_rate_hetero,'\n') ############################## # Courbes ROC & AUC ############################## seuil=seq(from=0.9, to=0.1,by=-0.01) # le seuil se déplace de 0,9 à 0,1 par centièmes tvp <- tfp <- NULL # Classifieur gaussien homoscédastique for (i in 1:length(seuil)){ fp=sum(testlabels=='healthy'&prediction_homo[6][,2]=seuil[i]) # nbre de faux négatifs vn=sum(testlabels=='healthy'&prediction_homo[6][,2]>=seuil[i]) # nbre de vrais négatifs tvp[i]=vp/(vp+fn) # taux de vrais positifs tfp[i]=fp/(fp+vn) # taux de faux positifs } plot(tfp,tvp,'l',col='blue', main='courbes ROC des classifieurs gaussiens',xlab='TFP',ylab='TVP') # courbe ROC du classifieur homoscédastique library(flux) AUC <- auc(tfp,tvp) cat('AUC pour ROC-classifieur homoscédastique : ',AUC,'\n') # Classifieur gaussien hétéroscédastique for (i in 1:length(seuil)){ fp=sum(testlabels=='healthy'&prediction_hetero[6][,2]=seuil[i]) # nbre de faux négatifs vn=sum(testlabels=='healthy'&prediction_hetero[6][,2]>=seuil[i]) # nbre de vrais négatifs tvp[i]=vp/(vp+fn) # taux de vrais positifs tfp[i]=fp/(fp+vn) # taux de faux positifs } points(tfp,tvp,'l',col='red') # courbe ROC du classifieur hétéroscédastique AUC <- auc(tfp,tvp) cat('AUC pour ROC-classifieur hétéroscédastique : ',AUC,'\n') legend('bottomright', c('homoscédastique','hétéroscédastique'), lty=c(1,1), col=c('blue','red') )