###############################################################################
## Inhalt:
## Kapitel 7.2 (Elemente Multivariater Statistik) des Kurses von 
## Ruckdeschel & Kohl
###############################################################################


###############################################################################
## Lesen Sie die Abschnitte 7.2.1 - 7.2.3 im Skript von Ruckdeschel & Kohl.
###############################################################################


###############################################################################
## Die multivariate Normalverteilung
###############################################################################
## Die folgenden Pakete enthalten Funktionen im Zusammenhang mit der
## multivariaten Normalverteilung:
## mvtnorm, mnormt, sn, fMultivar, monomvn, norm, splus2R, ...

###########################################################
## Die Dichte der multivariaten Normalverteilung
###########################################################
library(mvtnorm)
?dmvnorm

#######################################
## Plot der Dichte einer 2-dimensionalen Normalverteilung
#######################################
## 1. Setzen der Parameter
mittel <- c(1,2) # Mittelwertvektor
sigma <- matrix(c(4,2,2,3), ncol = 2) # Kovarianz

## 2. Vorgabe des Gitters und Berechnung der Dichte
x.grid <- seq(-5, 10, length = 100)
y.grid <- seq(-4, 8, length = 100)
xy.grid <- expand.grid(x.grid, y.grid)
dichte <- matrix(dmvnorm(xy.grid, mean = mittel, sigma = sigma), ncol = 100)

## 3. Farbauswahl für den Plot
library(RColorBrewer)
mypalette <- colorRampPalette(brewer.pal(9, "YlOrRd")) 

## 4. Plot mit Hilfe von image und contour
image(x = x.grid, y = y.grid, z = dichte, col = mypalette(20), xlab = "x",
      ylab = "y", main = "Dichte einer 2-dim. Normalverteilung")
contour(x = x.grid, y = y.grid, z = dichte, nlevels = 20, add = TRUE)

## 5. Plot mit Hilfe von persp
persp(x = x.grid, y = y.grid, z = dichte, col = brewer.pal(3, "Blues")[3], 
      phi = 30, theta = 30, xlab = "x", ylab = "y", 
      main = "Dichte einer 2-dim. Normalverteilung")


#######################################
## (Pseudo-)Zufallszahlen von einer 2-dim. Normalverteilung
#######################################
## 1. Setzen der Parameter
mittel <- c(1,2) # Mittelwertvektor
sigma <- matrix(c(4,2,2,3), ncol = 2) # Kovarianz

## 2. Simulation der (Pseudo-)Zufallszahlen
x <- rmvnorm(n = 10000, mean = mittel, sigma = sigma)

## 3. Standardplot
## aufgrund der großen Stichprobe "unübersichtlich"
plot(x, pch = 20, xlab = "x", ylab = "y", 
     main = "Bivariat-normalverteilte (Pseudo-)Zufallszahlen")


## 4. Plot inkl. Dichteschätzung der Punktewolke
## Mittels der Parameter "nrpoints" und "colramp" kann man steuern, ob und wie 
## die Datenpunkte eingezeichnet werden und welche Farbverläufe die Punktwolke 
## zeigen soll; z.B.
library(geneplotter)
smoothScatter(x, nrpoints = 10, 
              colramp = colorRampPalette(brewer.pal(9, "YlGnBu")), pch = 20)


#######################################
## Maximum Likelihood (ML) Schätzung im Fall der multivariaten Normalverteilung
#######################################
## 1. Setzen der Parameter (unbekannt, zu schätzen!)
mittel <- c(1,2) # Mittelwertvektor
sigma <- matrix(c(4,2,2,3), ncol = 2) # Kovarianz

## 2. Simulation der (Pseudo-)Zufallszahlen
x <- rmvnorm(n = 100, mean = mittel, sigma = sigma)

## 3. Schätzung für "mean" (Erwartungswert)
colMeans(x)

## 4. Schätzung für "sigma" (Kovarianz)
var(x)
## bzw.
cov(x)

## 5. Graphische Darstellung der Daten inkl. 95%-Konfidenzellipsoid
library(ellipse)
Kellipse <- ellipse(cov(x), centre = colMeans(x))
plot(Kellipse, type = "l", col = "orange", xlab = "x", ylab = "y", 
     main = "95%-Konfidenzellipsoid", lwd = 2)
points(x, pch = 20)


###############################################################################
## Hauptkomponenten- und Faktoranalyse
## Lesen Sie Abschnitt 7.2.4 im Skript von Ruckdeschel & Kohl
###############################################################################

###########################################################
## Hauptkomponnentenanalyse (PCA: principal component analysis)
###########################################################
## 1. Lade Datensatz
data(USArrests)

## 2. Berechne PCA 
pc11 <- princomp(USArrests) # basiert auf "eigen"
str(pc11) # Struktur des Ergebnisses
## vgl.
eigen(cov(USArrests))

## bzw.
pc12 <- prcomp(USArrests) # basiert auf "svd" - höhere Genauigkeit
str(pc12) # Struktur des Ergebnisses

## 3. Ausgabe der Ergebnisse
summary(pc11)
loadings(pc11) # die leeren Einträge sind nicht wirklich leer bzw. 0, sondern nur klein
## bzw.
pc12
summary(pc12)