#===========================================================================
# Zusatzmaterial zu Übungsblatt 2
# Statistik 4 für Nebenfachstudierende, SoSe 2017
#
# Datum: 04.05.2017
# Autoren: Nora Fenske , Margret Oelker, Valentin Wimmer, Micha Schneider
#=============================================================================


#=============================================================================
# Aufgabe 4: Multivariate Normalverteilung
#=============================================================================

# Gegeben:
X <- matrix(c(-.85, -2.8, -1.56, 1.62, -.85, .43), ncol=2, byrow=F)

#-----
# (a)
#-----
# Berechne Mittelwertsvektor
mean(X[,1])
mean(X[,2])

xquer2 <- c(mean(X[,1]), mean(X[,2]))
xquer2

# Alternative
xquer <- apply(X, MARGIN=2, FUN=mean)
xquer

#-----
# (b)
#-----
# Berechne zentrierte Datenmatrix

# manuell
Xstern <- cbind(X[,1] - xquer[1], X[,2] - xquer[2])
Xstern <- round(Xstern, digits=2)
Xstern


# Alternative: Definiere Matrix H
n <- 3
H <- diag(rep(1,n)) - 1/n*(rep(1,n)%*%t(rep(1,n)))
H

Xstern2 <- H%*%X
Xstern2 <- round(Xstern2, digits=2)
Xstern2

# Vergleiche beide Lösungen:
all(Xstern==Xstern2)


# Weitere Zentrierung ändert zentrierte Matrix nicht:
round(H%*%Xstern2, digits=2)
round(H%*%H%*%X, digits=2)


#-----
# (c)
#-----
# Berechne Kovarianzmatrix

# Ganz einfach mit R-Funktion
S <- cov(X)
S

# 1. Alternative: Berechne Werte einzeln
s11 <- 1/2 * sum((X[,1] - xquer[1])^2)
s22 <- 1/2 * sum((X[,2] - xquer[2])^2)
s12 <- 1/2 * sum((X[,1] - xquer[1])*(X[,2] - xquer[2]))

s11
s22
s12

# 2. Alternative: mit allgemeiner Matrixformel
x1 <- X[1,]
x2 <- X[2,]
x3 <- X[3,]

S2 <- 1/2 * ((x1 - xquer)%*%t(x1 - xquer) +
             (x2 - xquer)%*%t(x2 - xquer) +
             (x3 - xquer)%*%t(x3 - xquer))
S2
             
# 3. Alternative: mit Matrixformel für Varianzzerlegung
S3 <- 1/2 * ( x1%*%t(x1) + x2%*%t(x2) + x3%*%t(x3) -
              3 * xquer%*%t(xquer))
S3

#-----
# (d)
#-----
# Berechne Korrelationsmatrix

# Ganz einfach mit R-Funktion
R <- cor(X)
R

# Alternative: Wandle Kovarianz- in Korrelationsmatrix um
R2 <- cov2cor(S)
R2


#-----
# (e)
#-----
# Abbildungen der Dichte

# Paket für multivariate Normalverteilung laden
library(mvtnorm)

# Plot-Funktion
plotmvt <- function(covmat, muvec=c(0,0), xrange=c(-3,3), yrange=c(-3,3),
                    lgrid=50, perspPlot=TRUE){

    xgrid <- seq(xrange[1],xrange[2],length=lgrid)
    ygrid <- seq(yrange[1],yrange[2],length=lgrid)

    zmat <- matrix(0, nrow=lgrid, ncol=lgrid)
    for(i in 1:lgrid){
         for(j in 1:lgrid){
             zmat[i,j] <- dmvnorm(c(xgrid[i], ygrid[j]), mean=muvec, sigma=covmat)
         }
    }

    if(perspPlot==TRUE){
        par(mar=c(2,2,1,1))
        persp(x=xgrid, y=ygrid, z=zmat, xlab="", ylab="", zlab="",
              theta=-20,phi=30, cex.axis=1.15, ticktype="detailed")
    }
    else{
        par(mar=c(2,2,1,1))
        contour(x=xgrid, y=ygrid, z=zmat, cex.axis=1.15, labcex=1.2, xlab="", ylab="")
    } 

}

# Definiere wahren zugrunde liegenden Mittelwertsvektor und Kovarianzmatrix
meanvec <- c(-1,1)
covmat <- matrix(c(1, 0.9, 0.9, 2), ncol=2)

layout(matrix(c(1,2),ncol=2))
# Dichteplot
plotmvt(covmat=covmat, muvec=meanvec, lgrid=40, xrange=c(-3.5,1.5), yrange=c(-3,5))

# Konturplot
plotmvt(covmat=covmat, muvec=meanvec, lgrid=40, xrange=c(-3.5,1.5), yrange=c(-3,5),
        perspPlot=FALSE)


# falsche Korrelation (+ falscher Mittelwertsvektor)
meanvec <- c(1,1)
covmat <- matrix(c(5, -2, -2, 5), ncol=2)

plotmvt(covmat, lgrid=40, xrange=c(-5,6), yrange=c(-4,7), muvec=meanvec)
plotmvt(covmat, lgrid=40, muvec=meanvec, xrange=c(-5,6), yrange=c(-4,7), perspPlot=FALSE)

# falscher Mittelwert (aber wahre Kovarianzstruktur!)
meanvec <- c(-5,5)
covmat <- matrix(c(1, 0.9, 0.9, 2), ncol=2)

plotmvt(covmat, lgrid=40, xrange=c(-8,-1), yrange=c(0,10), muvec=meanvec)
plotmvt(covmat, lgrid=40, xrange=c(-8,-1), yrange=c(0,10),  muvec=meanvec, perspPlot=FALSE)
layout(matrix(1))

#-----
# Zusatz
#-----

# Ziehe 100 bzw. 1000000 Zufallszahlen aus der gleichen Dichte
X1 <- rmvnorm(100, mean=meanvec, sigma=covmat)
X2 <- rmvnorm(1000000, mean=meanvec, sigma=covmat)

# Berechne Mittelwertsvektor und Kovarianzmatrix
apply(X1, MARGIN=2, FUN=mean)
apply(X2, MARGIN=2, FUN=mean)
cov(X1)
cov(X2)

# Vergleiche mit der Wahrheit
meanvec
covmat


#=============================================================================
# Aufgabe 5: Zehnkampf
# Autoren: R Core Team, Valentin Wimmer, Micha Schneider
#=============================================================================

# auf der Diagonalen univariate Histogramme
panel.hist <- function(x, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5) )
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks; nB <- length(breaks)
  y <- h$counts; y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], y, col="lightgrey")
}

# im oberen Dreieck die Korrelationen als Text ausgeben
panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y,use="pairwise.complete.obs"))
  txt <- format(c(r, 0.123456789), digits=digits)[1]
  txt <- paste(prefix, txt, sep="")
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * (1+abs(r/3)))
}


mypoints <- function(x,y, ...) points(x,y, pch=".",cex=0.2)

# Daten laden
load("decathlon.Rdata")

# Graphik erstellen
pairs(decathlon[,8:17],
      diag.panel = panel.hist, cex.labels=1,upper.panel=panel.cor,
      panel=mypoints)

# Anmerkung: Graphik nicht identisch mit dem Aufgabenblatt, da dort 
#  eine zufällige Stichprobe aus den Originaldaten verwendet wurde.