rm(list=ls())

# 1/
# On télécharge le fichier olympic2.txt à partir de ma page web
#Il faut préciser au logiciel R, le répertoire qui contient ce fichier à l'aide de la commande setwd

setwd('C:/Users/christine/Desktop')
Olympic=read.table('olympic2.txt')
head(Olympic)		#cela permet de visualiser les 6 premières lignes

#on voit que l'importation des donnees est mauvaise car la première ligne contient le nom des variables

Olympic=read.table('olympic2.txt',header=TRUE)
Olympic=as.matrix(Olympic)		#afin de transformer le data.frame en une matrice

# on va faire l'analyse en composantes principales mais avant on va regarder un peu chacune des variables

for (i in 1:10)
{
X11()
boxplot(Olympic[,i])
}

#on remarque que certaines variables comportent des observations potentiellement abérantes.
#en regardant plus particulièrement la dernière observation, on s'aperçoit que cet athlète est mauvais
#dans toutes les épreuves et de loin
#donc on décide de le retirer de l'étude

Olympic=Olympic[-34,]

#par ailleurs, seules les variables relatives aux différentes épreuves nous intéressent
#donc on retire la colonne score

Olympic=Olympic[,-11]


# 2/
# a) dans toute ACP, le recentrage des données est obligatoire. Cela consiste à retirer à toutes les
#observations d'une variable, la moyenne de cette variable.

p=ncol(Olympic) # nb de variables
n=nrow(Olympic) # nb d'individus

Mcol=apply(Olympic,2,mean)	#vecteur ligne qui contient la moyenne de chacune des colonnes

#de façon équivalente
Mcol=c()
for (i in 1:p)
{
a=mean(Olympic[,i])
Mcol=c(Mcol,a)
}

#l'intérêt du apply est que l'on conserve le nom des colonnes dans le vecteur moyenne
#donc, c'est plus facile à comprendre!


# pour faire une grosse matrice avec la meme ligne repetee n fois
Mcolmat=matrix(rep(Mcol,n),ncol=p,byrow=TRUE)

#le rep(Mcol,n) permet de créer un vecteur où le vecteur Mcol est répété n fois
#l'option byrow=TRUE permet de remplir une matrice ligne par ligne et non pas colonne par colonne comme par défault


Oly_centre= Olympic-Mcolmat # Matrice des donnees recentrees

#pour vérifier que le recentrage est bon, il suffit de vérifier que la moyenne des colonnes de la nouvelle matrice
#est nulle pour toutes les colonnes

apply(Oly_centre,2,mean)
#on ne trouve pas 0 exactement mais toutes les valeurs sont très petites, de l'ordre du zéro du logiciel


# b) Ensuite on renormalise
# Attention : R calcule la version sans biais de la variance empirique (ou de l'ecart type empirique) en divisant par n-1 au lieu de n.
#par contre, la normalisation dans la fonction princomp se fait bien avec la version en 1/n

# Donc

SDcol=sqrt(n-1)/sqrt(n)*apply(Oly_centre,2,sd)	#pour avoir la version en 1/n
ones=as.matrix(rep(1,n)) #cela permet de créer un vecteur colonne ne contenant que des 1
SD_mat=ones%*%SDcol  #cela permet de créer une matrice où le vecteur SDcol est répété sur chacune des n lignes

Oly_renorm= Oly_centre/SD_mat # Matrice des donnees renormalisees

# 3/ maintenant on calcule l'acp
# a)
acp_olympic=princomp(Oly_renorm, scores=TRUE)

#cette commande est équivalente à demander de faire une ACP centrée réduite sur les données Olympic directement
#pour faire cela, il suffit de demander princomp(Olympic,cor=TRUE)
#le fait de mettre cor=TRUE oblige le logiciel à faire une ACP normée
#par défault, le logiciel fait une ACP non normée

#Rq de même princomp(Olympic) et princomp(Oly_centre) produit la même chose

# b)
(1/n)*(t(Oly_renorm)%*%Oly_renorm)
matcor=(n-1)/n*cov(Oly_renorm)

# Ces deux commandes renvoient la meme chose. C'est la matrice de covariance associee a Oly_renorm, avec la version 1/n.
#Rq: comme les données sont centrées réduites, la matrice de variance covariance de Oly_renorm est aussi la matrice de correlation
#Rq: c'est aussi la matrice de correlation des données Olympic
#Rq: la matrice de variance-covariance en 1/n de Olympic ((n-1)/n*cov(Olympic)) est 1/n*t(Oly_centre)%*%Oly_centre

# c)
summary(acp_olympic)

# Cette commande renvoie trois vecteurs. On a dans l'ordre:
# les racines carrees des valeurs propres de la matrice de correlation;
valpro=eigen(matcor)$values	#valeurs propres de la matrice de correlation de Olympic
sqrt(valpro)			#on retrouve la première ligne du summary

# la variance relative cad lambda_j / somme totale des lambda_k où les lambda_k sont les valeurs propres;
Propvariance=valpro/sum(valpro)

# et aussi la somme cumulee de ca
cumsum(Propvariance)

# Remarque : on voit que pour atteindre les 80% il faut aller jusqu'a la composante 5

# d)
acp_olympic$loadings
# Il s'agit des vecteurs propres de la matrice de correlation de Olympic

# Verification
eigen(matcor)$vectors

# Remarque : Attention, acp_olympics$loadings n'affiche que les coordonnees "franchement non nulles" pour pouvoir interpreter les axes
# pour  voir tous les non nuls vraiment on tape par exemple pour la composante 10

acp_olympic$loadings[,10]


# e)
acp_olympic$scores

# Renvoie les coordonnees des performances centrees reduites dans la nouvelle base (celle des vecteurs propres de la matrice de covariance)

# Verification :
Oly_renorm[1,]%*%acp_olympic$loadings	#ceci fournit la projection de des données centrées réduites de l'individu 1 sur les différentes composantes

# a comparer avec
acp_olympic$scores[1,]

#c'est la même chose

# f) Le diagramme des individus dans le premier plan factoriel de l'acp s'obtient de la maniere suivante

# Fait le plot  mais on sait pas ou sont les individus
plot(acp_olympic$scores[,1],acp_olympic$scores[,2],xlab="Comp.1",ylab="Comp.2") 
abline(h=0,v=0) # trace les deux axes au milieu pour qu'on y voit plus clair

# Comme propose par l'indication, si on veut voir les noms des individus (ici 1 2 ... )
plot(acp_olympic$scores[,1],acp_olympic$scores[,2],type='n',xlab="Comp.1",ylab="Comp.2") 
abline(h=0,v=0)
text(acp_olympic$scores[,1],acp_olympic$scores[,2],labels=1:33)

# Remarque : le premier axe est engendre principalement (cf la sortie de loadings) par le 110m haies et le 100m en positifs. 
# Les plus mauvais temps, ie les plus lents a ces deux epreuves sont bien les concurrents 17, 30 31 32 33 28 ... Bon c'est qualitatifs surtout que le premier axe a plein de composantes ...


# 4/
# a)

new100m=cor(Oly_renorm[,1],acp_olympic$scores)
# Sa norme euclidienne est donnee par
sqrt(sum(new100m^2))
# new100m est le vecteur des coordonnees de l'ancienne première variable "100m" (renormalisee) dans les nouvelles coordonnées

# Remarque : La norme euclidienne de new100m est egale a 1, donc il appartient a la sphere unite. 
# De plus, sa projection sur les deux premiers axes de l'ACP est de norme plus petite que 1 et appartient donc au disque de rayon 1

#pour comparer à la norme initiale de oly_renorm[,1], il faut calculer 1/sqrt(n)*sqrt(sum(oly_renorm[,1]^2)). Et oui, dans l'espace initial, les individus se voient chacun associer un poids 1/n

# b)
cor(Oly_renorm,acp_olympic$scores[,1])	#coordonnées des anciennes variables sur la première nouvelle composante
acp_olympic$sdev[1]*acp_olympic$loadings[,1]
# Ces deux vecteurs sont egaux. C'est une propriété de l'ACP
# Ils correspondent aux coordonnees des dix anciennes variables ("100m","long",...) sur le premier nouvel axe.

# c)
axe1=acp_olympic$sdev[1]*acp_olympic$loadings[,1] # coordonees sur axe 1 de toutes les anciennes variables
axe2=acp_olympic$sdev[2]*acp_olympic$loadings[,2] # coordonnees sur axe 2 de toutes les anciennes varaibles

plot(axe1,axe2,xlim=c(-1,1),ylim=c(-1,1),type='n',xlab="Comp.1 ",ylab="Comp.2")
t=seq(0,2*pi,0.01)	#calcul des coordonnées polaires de points du cercle unitaire
x=cos(t)
y=sin(t)
lines(x,y)		#ceci trace le cercle de corrélation qui est le cercle unitaire dans le cas des données centrées réduites
abline(h=0,v=0)
text(axe1,axe2,colnames(Olympic))

# Remarque : plus on est proche du cercle plus on est bien represente par ce qui se passe uniquement sur les nouvelles variables 1 et 2. 
# Typiquement ici hauteur etait peu important pour la matrice de covariance et se trouve releguee dans des composantes correspondants a des valeurs propres plus petite
# La premiere nouvelle variable ressemble surtout au 100 metre et 110 metres haies et a -perche et -longueur, 2 eme variable en gros "disque, 1500 m"


# 5
# a)
# les 4 plots

par(mfrow=c(2,2))
plot(acp_olympic$scores[,1],acp_olympic$scores[,2],type='n',xlab="Comp.1 ",ylab="Comp.2") 
abline(h=0,v=0)
text(acp_olympic$scores[,1],acp_olympic$scores[,2],labels=1:33)

plot(axe1,axe2,xlim=c(-1,1),ylim=c(-1,1),type='n',xlab="Comp.1 ",ylab="Comp.2")
t=seq(0,2*pi,0.01)
x=cos(t)
y=sin(t)
lines(x,y)
abline(h=0,v=0)
text(axe1,axe2,colnames(Olympic))

biplot(acp_olympic)

# le biplot ok fait tout en un, on peut retrouver les deux premiers superposes sur le troisieme 
# Attention : si on regarde bien les graduations rouges, on voit que les fleches tracees ne sont pas toutes dans le cercle unite
# Le cercle des correlations n'est pas de rayon 1 dans le biplot mais de rayon sqrt(n)

# b)
plot(acp_olympic,main="valeurs propres de l'ACP")
# Affiche les valeurs propres decroissantes de la matrice de covariance
# Peut servir a decider du choix du nombre d'axes a conserver

# c)
par(mfrow=c(1,2))
biplot(acp_olympic)
acp_brute=princomp(Olympic)	#cadre de l'acp non normée
biplot(acp_brute)

# Remarque : Sur les donnees brutes, les performances du 1500m (qui ont une variance beaucoup plus grande) ecrasent tout le reste
# Or, ce n'est pas le cas en vrai car le calcul des points du decathlon prend cela en compte

# d)
par(mfrow=c(1,2))
biplot(acp_olympic)
acp_easy=princomp(Olympic,cor=TRUE)
biplot(acp_easy)

# On voit que l'etape de Pre-traitement (Section 2) n'est pas necessaire
# On peut se contenter d'appliquer la fonction princomp avec l'option cor=TRUE pour que le logiciel applique la meme renormalisation


#2) Clustering
rm(list=ls())

#1/

data(iris)

iris

# iris est un data.frame , il contient pour 150 plantes, la longuer, largeur des sepales et des petales ainsi que l'espece d'iris correspondante

#2/
A=as.matrix(iris[,-5])

#3/

K=kmeans(A,3,iter.max=1,nstart=1)

# c'est la procedure kmeans, elle donne un clustering en 3 groupes et donne le centroide de chaque groupe, on associe un element au centroide le plus proche.

#4/

K$cluster #est un vecteur qui contient le numero du groupe pour chacune des
		#observations. Ainsi, l'element i est le numero du groupe de l'observation i.
		
K$centers # est une matrice a 3 lignes (on a demande 3 groupes). La ligne i contient les coordonnees du centre de gravite du groupe i.		

K$totss #contient la somme totale des carres des distances avec la moyenne globale a savoir sum_{pour i allant de 1 a n}t(x_i-xbar)*(x_i-xbar) 

#verification
xbar=colMeans(A)
a=0
n=nrow(A)
for(i in 1:n)
{
	a=a+ sum((A[i,]-xbar)^2)
}
abs(a-K$totss) #c'est la meme chose a la precision machine

K$withinss # vecteur des somme des carres des distances avec la moyenne du groupe pour chaque groupe

#verification
indice=which(K$cluster==1)
xcenter1=colMeans(A[indice,])
a=0
for(i in indice)
{
	a=a+ sum((A[i,]-xcenter1)^2)
}

abs(a-K$withinss[1]) #c'est la meme chose a la precision machine

K$tot.withinss # c'est la somme de tous les ecarts des points a leur centre de groupe respectif " somme des carres intragroupe"

# verification

K$tot.withinss==sum(K$withinss) # ca marche

K$betweenss # c'est la  " somme des carres intergroupe" cad totss-tot.withins

# verif


K$betweenss==(K$totss-K$tot.withinss)

K$size # nb de points dans chaque cluster

K$iter # nb d'iterations faite
K$ifault #pour expert de R, c'est marque dans l'aide. on ne s'en sert pas ...

#5/

# iter.max nb maximal d'iteration ici c'est 1, donc on ne fait qu'une iteration et ca n'a pas le temps de converger, d'ou le message d'erreur

# nstart nombre de reinitialisation au hasard et on prend le meilleur . ici 1 ce n'est pas bon car on va rester dans un minimum local

#6/

K=kmeans(A,3,iter.max=100,nstart=100)


# on peut controler sur le nombre d'iterations faites que cela converge tres vite.


#7/

clus1=which(K[[1]]==1)		#vecteur contenant le numero des observations du cluster 1
clus2=which(K[[1]]==2)		#vecteur contenant le numero des observations du cluster 2
clus3=which(K[[1]]==3)		#vecteur contenant le numero des observations du cluster 3

rep1=iris[clus1,5]			#valeur de la variable espece pour chacune des observations du cluster 1
rep2=iris[clus2,5]			#valeur de la variable espece pour chacune des observations du cluster 2
rep3=iris[clus3,5]			#valeur de la variable espece pour chacune des observations du cluster 3

	#on constate que si un des clusters ne contient que des observations associees a l'espece setosa
		#les deux autres clusters melangent des observations des especes versicolor et virginica. 
		
		#Comment expliquer la difference? 
	#une explication de cette difference est la suivante : si les 4 criteres morphologiques
		#permettent de bien differencier l'espece setosa des deux autres especes, elles ne 
		#permettent pas de distinguer clairement les especes versicolor et virginica. En effet
		#une plante n'est pas uniquement caracterisee  par ces 4 criteres.

#8/

n=nrow(A)-1 # kmeans veut toujours strictement plus de donnees que de groupes 
Kvec=1:n
v=rep(0,n)
for(k in Kvec)
{
	K=kmeans(A,k,iter.max=100,nstart=100)
	v[k]=K$tot.withinss
}

plot(Kvec,v)

# on retrouve bien que plus on va faire un grand nombre de groupes et plus la 
#variance intra-groupe va etre faible. Cependant, s'il y a une forte decroissance
		#au depart, ensuite la decroissance est plus lente. Il faut donc faire un
		#compromis entre une faible variance intra-groupes et un nombre de groupe raisonnable. Ici ca a l'air d'etre deux ou 3 groupes
	
	# 9 /
	
D=dist(A)
Db=dist(A,method="maximum")

# D est une "sorte" de matrice trinagulaire inferieure de class "dist" pour en faire une matrice symmetrique avec diagonale de 0 on fait as.matrix(D). Dans D sont calculees les distances euclidiennes entre paires de points et dans Db la meme chose avec la distance en norme infinie

#10/

hc=hclust(D^2,method='ward.D') # fait tout avec la distance D2
plot(hc)

# cela fait une classification hierarchique basee sur la matrice des distance euclidienne couple a la methode de ward	. Puis ca trace le dendogramme associe

hc$merge	# cf l'aide ca dit comment les paquets fusionnent dans le dendogramme au fur et a mesure -j "nouvelle entree d'une observation +k paquet cree a l'iteration k

hc$height # cf l'aide c'est ce qui est trace en hauteur dans le dendogramme. Pour la methode de ward c'est l'augmentation d'inertie intragroupe par le fait d'avoir fusionne. Dans le cas de Ward c'est la somme des carres intragroupe renormalise par le nombre total

# verification		
		
Dmat=as.matrix(D)		
		
Dmat[102,143]
Dmat[8,40] # premeire paire non nulle et dans 
hc$height[2] # j'ai bien le carre de la premiere paire non egale mergee

# a l'etape 22 on merge 41 et (1,18) ie deuxieme paire

# la formule de cours explique tout http://pbil.univ-lyon1.fr/R/pdf/stage7.pdf (p35)

# ici ca donne 

2/3*Dmat[1,41]^2+2/3*Dmat[18,41]^2-1/3*Dmat[1,18]^2 

hc$height[23] 

# si onfait

hcnew=hclust(D,method='ward.D2') 
plot(hcnew)

# on a alors que

hc$merge==hcnew$merge # meme arbre

abs(sqrt(hc$height)-hcnew$height)
# les longueurs sont maintenant vraiment des distances et pas les distances au carre

# dans les deux cas clairement plus grand saut a deux composantes ... il faudrait prendre k=2 pour les kmeans

K2=kmeans(A,2,iter.max=100,nstart=100)

clus1=which(K2[[1]]==1)		
clus2=which(K2[[1]]==2)

rep1=iris[clus1,5]			
rep2=iris[clus2,5]	

#on ne retrouve pas parfaitement setosa

# on pourrait prendre  l'acp et colorie pour voir si on y voit quelque chose ...