# Simulationen: Vorlesungsbeispiel

#-----------------------------------------------
# Daten generieren
# einen normalverteilten Zufallsvektor ziehen, mit n_obs Beobachtungen

# Parameter
# mu:     wahrer Mittelwert
# sigma:  wahre Standardabweichung
# n_obs:  Anzahl an Beobachtungen

simulate_data <- function(mu, sigma, n_obs){
  x <- rnorm(n=n_obs, mean=mu, sd=sigma)
  x
}

# Beispielaufruf
testdata <- simulate_data(mu=50, sigma=10, n_obs=100)


#-----------------------------------------------
# generierte Daten analysieren
# den Mittelwert der simulierten Beobachtungen berechnen

# Parameter
# dat: Datenvektor

analyze_data <- function(dat){
  mean_x <- mean(dat)  
  mean_x
}

# Beispielaufruf
mean <- analyze_data(dat=testdata)


#-----------------------------------------------
# Beurteilungskriterium berechnen 
# Absolutbetrag der Differenz von wahrem Mittelwert und geschätzem Mittelwert

# Parameter
# mu_est: geschätzter Mittelwert
# mu:     wahrer Mittelwert 

analyze_result <- function(mu_est, mu){
  diff <- abs(mu-mu_est)
  diff 
}

# Beispielaufruf
mean_diff <- analyze_result(mu_est=mean, mu=50)


#-----------------------------------------------
# beides verknüpfen und mehrmals aufrufen

# Parameter
# mu:             wahrer Mittelwert
# sigma:          wahre Standardabweichung
# n_obs:          Anzahl an Beobachtungen
# n_repetitions:  Anzahl an Wiederholungen

# ohne seed (aber seed für Gesamtsimulation sinnvoll)
one_simulation <- function(mu, sigma, n_obs, n_repetitions){
  
  # Ergebnisvektor initialisieren
  mean_diff <- rep(NA, times=n_repetitions)
  
  # Daten generieren, analysieren, Beurteilungskriterum berechnen und Ergebnisvektor füllen
  for(i in 1:n_repetitions){
  dat <- simulate_data(mu=mu, sigma=sigma, n_obs=n_obs)
  res <- analyze_data(dat=dat)
  diff <- analyze_result(mu_est=res, mu=mu)
  
  mean_diff[i] <- diff
  }
  
  # Rückgabewert
  mean_diff
}

# mit seed (zur Reproduzierbarkeit der Einzelergebnisse)
one_simulation <- function(mu, sigma, n_obs, n_repetitions, seed=1:n_repetitions){
  
  # Ergebnisvektor initialisieren
  mean_diff <- rep(NA, times=n_repetitions)
  
  # Daten generieren, analysieren, Beurteilungskriterum berechnen und Ergebnisvektor füllen
  for(i in 1:n_repetitions){
    set.seed(seed[i])
    dat <- simulate_data(mu=mu, sigma=sigma, n_obs=n_obs)
    res <- analyze_data(dat=dat)
    diff <- analyze_result(mu_est=res, mu=mu)
    
    mean_diff[i] <- diff
  }
  
  # Rückgabewert
  mean_diff
}


# Beispielaufrufe
one_simulation(mu=50, sigma=10, n_obs=100, n_repetitions=100)
one_simulation(mu=50, sigma=20, n_obs=100, n_repetitions=100)
one_simulation(mu=50, sigma=30, n_obs=100, n_repetitions=100)



#-----------------------------------------------
# mit variierenden Parametern aufrufen

# Parameter
# mu_vector:      Vektor mit wahren Mittelwerten
# sigma_vector:   Vektor mit wahren Standardabweichungen
# n_obs_vector:   Vektor mit variierenden Anzahlen an Beobachtungen
# n_repetitions:  Anzahl an Wiederholungen


complete_simulation <- function(mu_vector, sigma_vector, n_obs_vector, n_repetitions){
 
  # alle Möglichen Parameterkombinationen bilden
  parameter_combinations <- expand.grid(mu_vector,
                                        sigma_vector,
                                        n_obs_vector) 
  
  # sinnvolle Variablennamen definieren
  names(parameter_combinations) <- c("mu", "sigma", "n_obs") 
  # print(parameter_combinations)
  
  # Ergebnismatrix initialisieren
  result <- matrix(NA, ncol=nrow(parameter_combinations), nrow=n_repetitions)
  
  # Ergebnismatrix befüllen
  for(i in 1:ncol(result)){
    # cat("Iteration ist ", i, "\n")
    result[,i] <- one_simulation(mu=parameter_combinations$mu[i], 
                                 sigma=parameter_combinations$sigma[i], 
                                 n_obs=parameter_combinations$n_obs[i],
                                 n_repetitions=n_repetitions)
  }
  
  # Rückgabewert erstellen
  ret <- list(parameter_combinations=parameter_combinations,
              result=result)
  # Rückgabewert
  ret
}
  

#--------------------
# Funktionsaufruf (jetzt mit Startwert für Reproduzierbarkeit)
set.seed(8372091)
set.seed(9478205)
final <- complete_simulation(mu_vector=50, 
                             sigma_vector=c(10,20,30), 
                             n_obs_vector=c(10,100,1000), 
                             n_repetitions=100)
# Ergebnis speichern!
# save(final, ...)


#---------------------
# Auswertung
# Ergebnis anschauen und graphisch darstellen
final$parameter_combinations
boxplot(final$result)


#--------------------
# Zeitabschätzung (mit kleinerem Umfang)

# mit system.time 
# einmalige Ausführung und Rückgabe der Zeit in Sekunden
system.time(
complete_simulation(mu_vector=50, 
                    sigma_vector=c(10,20,30), 
                    n_obs_vector=c(10,100,1000), 
                    n_repetitions=50)
)

# mit package 'microbenchmark' 
# mehrmalige Ausführung mit summary über die benötigte Zeit, Zeiteinheit flexibel 
library("microbenchmark")
mb <- microbenchmark(complete_simulation(mu_vector=50, 
                                         sigma_vector=c(10,20,30), 
                                         n_obs_vector=c(10,100,1000), 
                                         n_repetitions=50), 
                     times=100, unit="s")
mb


#-----------------------------------------------
# Möglichkeiten, um einen Startwert zu wählen

# geeignete Zufallszahl ziehen und dann fixieren
round(runif(n=1, min=1000, max=100000))
fixed_seed <- 86889
set.seed(fixed_seed)
# ...


# Datum von heute
fixed_seed <- 11062014
set.seed(fixed_seed)
# ...