### Multivariate Verfahren Tutorium 1 

### Aufgabe 1

# n=3 Zeilen und p=2 Spalten
X <- matrix(nrow=3, ncol=2)
X[,1] <- c(3,1,2) # Spalte 1
X[,2] <- c(2,5,5) # Spalte 2

### a
# Spaltenmittelwerte
mean(X[,1])
mean(X[,2])

## schnell
# MARGIN=2 => wende Funktion mean() auf die Spalten von X an!
x_quer <- apply(X , 2, mean)
x_quer

### b (iii)
## mit Formel fuer erwartungstreuen Schaetzer aus Folie 15
n <- nrow(X)
sum_xi_txi <- X[1,]%*%t(X[1,]) + X[2,]%*%t(X[2,]) + X[3,]%*%t(X[3,])
S <- 1/(n-1) * (sum_xi_txi - n*x_quer%*%t(x_quer))
S

## mit Zentrierungsmatrix H aus Folie 16
H <- diag(n)-(1/n)*matrix(1,nrow=n,ncol=n)
S <- 1/(n-1)*t(X)%*%H%*%X
S

## schnell
cov(X)


### Aufgabe 2

set.seed(123)

#Definition der Kovarianzmatrix 
S<-matrix(c(2,2,2,5),nrow=2,ncol=2)

# zur Generierung bivariater normalverteilter Zufallszahlen
library(mvtnorm)

# zunaechst einfach: Generierung einer Wishart-verteilten Zufallsmatrix mit 
# Parameter S und m=1
x <- rmvnorm(1,c(0,0),S) # eine 2-dim normalverteile ZV 
M <- x[1,]%*%t(x[1,])    # eine Wishart-verteile ZV mit m=1 
M

#Funktion zur Berechnung von n Wishart-verteilten Zufallsmatrizen mit Parametern
#S und m
m<-10 # Freiheitsgrade
n<-100

rwishart<-function(n,S,m){ 
  
  mu1 <- rep(0,length(S[,1]))  #Definition von \mu=0
  #Dimension von Sigma liefert, dass eine 2-dimensionale
  #Multivariate Normalverteilung zugrunde liegt!
  
  #Simulation der normalverteilten Zufallsvektoren
  A<-rmvnorm(n*m,mu1,S) # simuliere 2-dim. normalvertielte Zufallszahlen mit Erwartungswert 0 und
  # Kovarianzmatrix Sigma
  # Benoetigt werden n*m Zufallszahlen, da fuer jede Wishart-verteilte
  # Zufallszahl m normalverteilte Zufallszahlen benoetigt werden!
  # dim(A)
  wlist<-list()  # Liste, weil die Zufallszahlen Matrizen sind (der Dimension = dim(S))
  #Und nun jeweils m xx' aufsummieren zu einer Wishart Zufallsvariable
  for(j in 1:n){
    W <- 0
    # fuer die erste wishart-verteilte Zufallszahl die ersten m normalverteilten Zufallszahlen
    # fuer die zweite wishart-verteilte Zufallszahl die zweiten m normalverteilten Zufallszahlen
    # ...
    # fuer die n-te wishart-verteilte Zufallszahl die letzten m normalverteilten Zufallszahlen
    for(i in ((j-1)*m+1):((j-1)*m+m)){
      W<-W+A[i,]%*%t(A[i,])
    }
    #Liste, die die Wishart-verteilten Zufallsvariablen enthaelt
    wlist[[j]]<-W
  }
  return(wlist)
}

# Ziehe 100 wishart-verteilte Zufallszahlen mit Kovarianzmatrix S und 10 Freiheitsgraden
wishartsample <- rwishart(n=100,S=S,m=10)
wishartsample[[1]]
wishartsample[[100]]
#Das ist eine Liste von W(S,m)-verteilten Matrizen.