# Stochastik IV, Übungsblatt 4

# Aufgabe 3

library(MASS)
library(ellipse)

# Teilaufgabe (a)

mu <- c(1, 2.5)
sigma <- matrix(c(2, -3, -3, 5), 2, 2)
cov2cor(sigma)  # Variablen sind sehr stark negativ korreliert

plot(ellipse(sigma, centre = mu, level = 0.99), type = "l", 
     main = "Contour Plot einer 2D Normalverteilung")
lines(ellipse(sigma, centre = mu, level = 0.95), col = "2")
lines(ellipse(sigma, centre = mu, level = 0.80), col = "3")
lines(ellipse(sigma, centre = mu, level = 0.50), col = "4")

points(1.0, 2.5, pch = 19)
text(1.25, 2.15, "(1, 2.5)")
legend(2, 8, c("0.99", "0.95", "0.80", "0.50"), col = c(1,2,3,4), lty = 1)

# Teilaufgabe (b)

windows()

# Einlesen der Daten und Einbinden in Suchpfad

data.olives <- 
read.table("D:/Univ Augsburg/Lehre/SS07/VO_MultiVar/Datensaetze/Olives.txt",
head = TRUE, sep = "\t", quote = "")

attach(data.olives)

# Geschätzte Größen

(mean.vec <- c(mean(palmitoleic), mean(oleic)))
(cov.matrix <- cov(cbind(palmitoleic, oleic)))
cov2cor(cov.matrix)  # Variablen sind stark negativ korreliert

# Contour Plot 

# Ellipsen der 2D Normalverteilungsdichte

plot(ellipse(cov.matrix, centre = mean.vec, level = 0.99), col = 2, lwd = 2,
     type = "l", 
     main = "Contour Plot von Kern- und Normalverteilungsdichteschätzer")
lines(ellipse(cov.matrix, centre = mean.vec, level = 0.95), col = 3, lwd = 2)
lines(ellipse(cov.matrix, centre = mean.vec, level = 0.80), col = 4, lwd = 2)
lines(ellipse(cov.matrix, centre = mean.vec, level = 0.50), col = 5, lwd = 2)
legend(160, 8500, c("Kerndichte", "Normaldichte 0.99", "Normaldichte 0.95", 
                    "Normaldichte 0.80", "Normaldichte 0.50"), 
       col = c(1,2,3,4,5), lty = 1)

# Höhenlinien des Kerndichteschätzers

points(palmitoleic, oleic, col = "lightblue")
dens.est.bandwidth.nrd <- kde2d(palmitoleic, oleic, 
                                h = c(bandwidth.nrd(palmitoleic), 
                                      bandwidth.nrd(oleic)), n = 30)
contour(dens.est.bandwidth.nrd, nlevels = 5, method = "flattest", lwd = 2,
        add = TRUE)