library(mvtnorm)

Simus <- 5000

# Total MSE! 
# (wichtig, den Gesamt-MSE zu betrachten, nicht komponentenweise )
MSE.TE <- vector(length=Simus, mode="numeric")
MSE.TEstar <- vector(length=Simus, mode="numeric")
MSE.TEtilde <- vector(length=Simus, mode="numeric")

for ( s in 1:Simus){

# Daten erzeugen (Varianz 1, ansonsten muss Formel veraendert werden!)
n <- 1
m <- 3
mu.true <- c(1,2,0.9)
x <- rmvnorm(n=n, mean=mu.true, sigma=diag(x=1,nrow=3) )


t <- colMeans(x)
MSE.TE[s] <- crossprod( t - mu.true )[1]
MSE.TEstar[s] <- crossprod((1- ((m-2)/n)/crossprod(t)[1] ) * t - mu.true)[1]
MSE.TEtilde[s] <- crossprod(max(0, (1- ((m-2)/n)/crossprod(t)[1] )) * t - mu.true)[1]
}

cat("Empirischer MSE des ?blichen Sch?tzers (Mittelwert): ", mean(MSE.TE), "\n")
cat("Empirischer MSE des Stein-Sch?tzers:                 ", mean(MSE.TEstar), "\n")
cat("Empirischer MSE des modifizierten Stein-Sch?tzers:   ", mean(MSE.TEtilde), "\n")


# Componentwise MSE
MSE.TE <- matrix(NA, ncol = m, nrow=Simus)
MSE.TEstar <- matrix(NA, ncol = m, nrow=Simus)
MSE.TEtilde <- matrix(NA, ncol = m, nrow=Simus)

for ( s in 1:Simus){
  
  # Daten erzeugen (Varianz 1, ansonsten muss Formel veraendert werden!)
  n <- 1
  m <- 3
  mu.true <- c(1,2,0.9)
  x <- rmvnorm(n=n, mean=mu.true, sigma=diag(x=1,nrow=3) )
  
  
  t <- colMeans(x)
  MSE.TE[s,] <- ( t - mu.true )^2
  MSE.TEstar[s,] <- ((1- ((m-2)/n)/crossprod(t)[1] ) * t - mu.true)^2
  MSE.TEtilde[s,] <- (max(0, (1- ((m-2)/n)/crossprod(t)[1] )) * t - mu.true)^2
}

cat("Empirische MSEs des ?blichen Sch?tzers (Mittelwert): ", colMeans(MSE.TE), "\n")
cat("Empirische MSEs des Stein-Sch?tzers:                 ", colMeans(MSE.TEstar), "\n")
cat("Empirische MSEs des modifizierten Stein-Sch?tzers:   ", colMeans(MSE.TEtilde), "\n")

cat("Empirischer MSE des ?blichen Sch?tzers (Mittelwert): ", sum(colMeans(MSE.TE)), "\n")
cat("Empirischer MSE des Stein-Sch?tzers:                 ", sum(colMeans(MSE.TEstar)), "\n")
cat("Empirischer MSE des modifizierten Stein-Sch?tzers:   ", sum(colMeans(MSE.TEtilde)), "\n")