###### Aufgabe 2 #########

# a)

#Einlesen des Datensatzes
euro <- read.table("https://www.elab.moodle.elearning.lmu.de/pluginfile.php/36839/mod_resource/content/1/europa.txt",header=T)
Xeuro <- data.frame(euro[,-1],row.names=euro[,1])

#standardisieren der Daten
euroscaled <- scale(Xeuro,scale=T,center=T)

# b) 

# quadrierte euklidische Distanzen
dist_eucl2 <- dist(euroscaled, method="euclidean")^2

# Single Linkage mit allen Kovariablen 
sl   <- hclust(dist_eucl2, method="single")

# c) 

# Zentroid mit allen Kovariablen 
cent <- hclust(dist_eucl2, method="centroid")

# d) 

# Complete Linkage mit allen Kovariablen und Mahalanobis-Distanz 
S <- var(euroscaled)
dist_mahal <- matrix(NA,nrow=24,ncol=24)
for(i in 1:24){
  dist_mahal[i,] <- mahalanobis(euroscaled, euroscaled[i,],S)
}
rownames(dist_mahal) <- colnames(dist_mahal) <- labels(dist_eucl2)
dist_mahal <- as.dist(dist_mahal)

cl <- hclust(dist_mahal, method="complete")

# e) 
par(mfrow=c(1,3))
plot(sl)
plot(cent)
plot(cl)


# f) 

#Das package cluster enthält verschiedene Clustering-Verfahren (sowohl hierarchische
#als auch Partitionsverfahren)
library(cluster)
help(kmeans)
#Partitionsverfahren, es muss Datenmatrix (standardisiert) und Anzahl der Cluster
#?bergeben werden (hier k=4. Ein erster Anhaltspunkt zur Clusteranzahl k?nnen hierarchische
#Verfahren liefern.)

# kmeans-Clustering mit allen Kovariablen, zehn zufällige Startpartitionen
set.seed(77)
km4 <- kmeans(x=euroscaled,centers=4,iter.max=20,nstart=10)

#Damit erhält man folgende Cluster:
km4.cl <- cbind(as.character(euro[order(km4$cluster),1]),sort(km4$cluster))


# g) 

#Zur Visualisierung beschränken wir uns auf die Variablen arbl und brut

# Single Linkage mit arbl und brut und k=4
sl2     <- hclust(dist(euroscaled[,3:4]),method="single")
sl24    <- cutree(sl2,k=4)
sl24.cl <- cbind(as.character(euro[order(sl24),1]),sort(sl24))

# Zentroid mit arbl und brut und k=4
cent2     <- hclust(dist(euroscaled[,3:4]),method="centroid")
cent24    <- cutree(cent2,k=4)
cent24.cl <- cbind(as.character(euro[order(cent24),1]),sort(cent24))

# Complete Linkage mit arbl und brut und k=4 und Mahalanobis-Distanz 
S <- var(euroscaled[,3:4])
dist_mahal <- matrix(NA,nrow=24,ncol=24)
for(i in 1:24){
  dist_mahal[i,] <- mahalanobis(euroscaled[,3:4], euroscaled[i,3:4],S)
}
dist_mahal <- as.dist(dist_mahal)
cl2     <- hclust(dist_mahal, method="complete")
cl24    <- cutree(cl2,k=4)
cl24.cl <- cbind(as.character(euro[order(cl24),1]),sort(cl24))

# kmeans-Clustering mit arbl und brut , zehn zufällige Startpartitionen
set.seed(77)
km2 <- kmeans(x=euroscaled[,3:4],centers=4,iter.max=20,nstart=10)
km2.cl <- cbind(as.character(euro[order(km2$cluster),1]),sort(km2$cluster))

# Graphischer Vergleich
plot_clust <- function(res,main=""){
  plot(euro[res==1,4],euro[res==1,5],type="n",xlim=c(0,37000),ylim=c(0,20),xlab="brut",ylab="arbl",main=main)
  for(i in 1:4){
    text(euro[res==i,4],euro[res==i,5],labels=as.character(euro[res==i,1]),col=i)
  }
}



par(mfrow=c(2,2))

# kmeans
plot_clust(km2$cluster,main="kmeans")

# single
plot_clust(sl24,main="single")

# centroid
plot_clust(cent24,main="centroid")

# complete
plot_clust(cl24,main="complete")

par(mfrow=c(1,1))



######## Aufgabe 3 ########


# b) 

library(MASS)
data(geyser)

## kurzer Plot:
plot(geyser)
# Man erkennt drei Klassen. Beachte den Rundungseffekt (viele Beob. mit duration = 2 oder 4)


# c) 

## Package zum Modellbasierten Clustern
library(mclust)

## Dokumentation:
help(Mclust)
help(mclustModelNames)


set.seed(0)
### erstes Modell: Cluster mit zum Koordinatensystem parallelen Achsen, gleicher Groesse:
m1 <- Mclust(data=geyser, G=3, modelNames = rep("EII", 3))
summary(m1, parameters=TRUE)                                 ## BIC = -3770.529


plot(m1, what="classification")        ## Daten mit berechneter Clusterung und Normalverteilungsellipse (Achsenl?nge = Standardabweichung)
plot(m1, what="uncertainty")           ## Unsicherheit: je dicker der Punkt, desto unsicherer ist die Zuteilung der Beobachtung zu einem Cluster
plot(m1, what="density")               ## Konturlinien der Mischverteilung    (vogelwild...)

## Zwischenfazit: die parallelen und gleich grossen Cluster sind hier ungeeignet


### zweites Modell: parallele Cluster mit unterschiedlicher Groesse und Form:
m2 <- Mclust(data=geyser, G=3, modelNames = rep("VVI", 3))
summary(m2, parameters=TRUE)                                 ## BIC = -2813.52

plot(m2, what="classification")       
plot(m2, what="uncertainty")           
plot(m2, what="density")                                     ## viel besser im Vergleich zu m1


### drittes Modell: nicht-parallele Cluster mit unterschiedlicher Groesse und Form:    (maximale Flexibilit?t)
m3 <- Mclust(data=geyser, G=3, modelNames = rep("VVV", 3))
summary(m3, parameters=TRUE)                                 ## BIC = -3059.029

plot(m3, what="classification")                              ## passt nicht
plot(m3, what="uncertainty")           
plot(m3, what="density")                                     ## passt nicht



#### m2 liefert mit Abstand das beste BIC (aller moeglicher Spezifikationen, auch der nicht Gezeigten).
#### Optisch scheint m2 ebenfalls sinnvoll zu clustern.

## geschaetzte posteriori-Wkt. aller Beobachtungen:
round(m2$z, 3)
## zu Diagnosezwecken betrachtet man haeufig folgendes Histogram, bei dem sich bei 
## gut separierten Mischkomponenten moeglichst alle Werte nahe 0 oder nahe 1 befinden sollten:
hist(round(m2$z, 3))

## Klassifikation:
m2$classification
## Unsicherheit bezueglich der Zuteilung:
round(m2$uncertainty, 3)
table(round(m2$uncertainty, 3))