#Aufgabe 2

# a)

#Einlesen des Datensatzes
euro <- read.table("http://www.statistik.lmu.de/institut/lehrstuhl/semsto/Lehre/Multivariate2014/europa.txt",header=T)
Xeuro <- data.frame(euro[,-1],row.names=euro[,1])

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

#Das package cluster enthält verschiedene Clustering-Verfahren (sowohl hierarchische
#als auch Partitionsverfahren)
library(cluster)


# b)

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))

#Zur Visualisierung beschränken wir uns auf die Variablen arbl und brut
set.seed(77)
km2 <- kmeans(x=euroscaled[,3:4],centers=4,iter.max=20,nstart=10)
km2

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

#Graphisch dargestellt:
par(mfrow=c(1,1))
plot(euro[km2$cluster==1,4],euro[km2$cluster==1,5],type="n",xlim=c(0,37000),ylim=c(0,20),xlab="brut",ylab="arbl")
text(euro[km2$cluster==1,4],euro[km2$cluster==1,5],labels=as.character(euro[km2$cluster==1,1]),col=1)
text(euro[km2$cluster==2,4],euro[km2$cluster==2,5],labels=as.character(euro[km2$cluster==2,1]),col=2)
text(euro[km2$cluster==3,4],euro[km2$cluster==3,5],labels=as.character(euro[km2$cluster==3,1]),col=3)
text(euro[km2$cluster==4,4],euro[km2$cluster==4,5],labels=as.character(euro[km2$cluster==4,1]),col=4)


# c) 

# Single Linkage mit allen Kovariablen und k=4
# sl4 <- hclust(dist(euroscaled),method="single")
# sl44 <- cutree(sl4,k=4)
# sl44.cl <- cbind(as.character(euro[order(sl44),1]),sort(sl44))

# 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 allen Kovariablen und k=4
# cent4 <- hclust(dist(euroscaled),method="centroid")
# cent44 <- cutree(cent4,k=4)
# cent44.cl <- cbind(as.character(euro[order(cent44),1]),sort(cent44))

# 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))


# Graphischer Vergleich
par(mfrow=c(1,3))

# kmeans
plot(euro[km2$cluster==1,4],euro[km2$cluster==1,5],type="n",xlim=c(0,37000),ylim=c(0,20),xlab="brut",ylab="arbl")
text(euro[km2$cluster==1,4],euro[km2$cluster==1,5],labels=as.character(euro[km2$cluster==1,1]),col=1)
text(euro[km2$cluster==2,4],euro[km2$cluster==2,5],labels=as.character(euro[km2$cluster==2,1]),col=2)
text(euro[km2$cluster==3,4],euro[km2$cluster==3,5],labels=as.character(euro[km2$cluster==3,1]),col=3)
text(euro[km2$cluster==4,4],euro[km2$cluster==4,5],labels=as.character(euro[km2$cluster==4,1]),col=4)

# single
plot(euro[sl24==1,4],euro[sl24==1,5],type="n",ylim=c(0,20),xlim=c(0,37000),xlab="brut",ylab="arbl")
text(euro[sl24==1,4],euro[sl24==1,5],labels=as.character(euro[sl24==1,1]),col=1)
text(euro[sl24==2,4],euro[sl24==2,5],labels=as.character(euro[sl24==2,1]),col=2)
text(euro[sl24==3,4],euro[sl24==3,5],labels=as.character(euro[sl24==3,1]),col=3)
text(euro[sl24==4,4],euro[sl24==4,5],labels=as.character(euro[sl24==4,1]),col=4)

# centroid
plot(euro[cent24==1,4],euro[cent24==1,5],type="n",ylim=c(0,20),xlim=c(0,37000),xlab="brut",ylab="arbl")
text(euro[cent24==1,4],euro[cent24==1,5],labels=as.character(euro[cent24==1,1]),col=1)
text(euro[cent24==2,4],euro[cent24==2,5],labels=as.character(euro[cent24==2,1]),col=2)
text(euro[cent24==3,4],euro[cent24==3,5],labels=as.character(euro[cent24==3,1]),col=3)
text(euro[cent24==4,4],euro[cent24==4,5],labels=as.character(euro[cent24==4,1]),col=4)

par(mfrow=c(1,1))






######## Aufgabe 3 ########
library(MASS)
data(geyser)

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

## Package zum Modellbasierten Clustern
library(mclust)

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


set.seed(0)
### erstes Modell: Cluster mit zum Koordinatensystem parallelen Achsen, gleicher Größe:
m1 <- Mclust(data=geyser, G=2, modelNames = rep("EII", 2))
summary(m1, parameters=TRUE)                                 ## BIC = -3917.92


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 großen Cluster sind hier ungeeignet


### zweites Modell: parallele Cluster mit unterschiedlicher Größe und Form:
m2 <- Mclust(data=geyser, G=2, modelNames = rep("VVI", 2))
summary(m2, parameters=TRUE)                                 ## BIC = -2897.027

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 Größe und Form:    (maximale Flexibilität)
m3 <- Mclust(data=geyser, G=2, modelNames = rep("VVV", 2))
summary(m3, parameters=TRUE)                                 ## BIC = -3030.993

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



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

## geschätzte posteriori-Wkt. aller Beobachtungen:
round(m2$z, 3)
## zu Diagnosezwecken betrachtet man häufig folgendes Histogram, bei dem sich bei 
## gut separierten Mischkomponenten möglichst alle Werte nahe 0 oder nahe 1 befinden sollten:
hist(round(m2$z, 3))

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