Skip to content
Snippets Groups Projects
Commit fda26fb1 authored by Jaroslav Borodavka's avatar Jaroslav Borodavka
Browse files

Translated to German and added a preamble.

parent a1abc750
No related branches found
No related tags found
No related merge requests found
#################################################################################################
#################################################################################################
##
#################################################################################################
#################################################################################################
# required packages
library(Directional)
library(MASS)
library(MonteCarlo)
library(data.table)
setwd("/Users/faik/Desktop/Mathematik/Masterarbeit/Ansatz_2/R-Codes")
setwd("/home/jaroslav-borodavka/Schreibtisch/Dissertation/Arbeit/Bernoulli Submission/codes/codes_git_project/R")
#################################################################################################
## Funktionen ##
## functions ##
# Kovarianzkern
# covariance kernel formulas as per section 2 of the main manuscript
rho_beta <- function(b, c, beta){
d = dim(b)[2]
res = switch (beta,
......@@ -22,51 +29,70 @@ rho_beta <- function(b, c, beta){
return(res)
}
# uniformly distributed vectors on S^(d-1)
runif_sphere <- function(n, d){
X = mvrnorm(n, rep(0,d), diag(1,d))
if (n == 1){
temp = as.vector(crossprod(X))
U = X/sqrt(temp)
}
else{
temp = diag(X %*% t(X))
U = matrix(0, nrow = n, ncol = d)
for (j in 1:n){
U[j,] = X[j,]/sqrt(temp[j])
}
}
return(U)
}
T_asymptotic_stat_value <- function(d, beta){
m = 0
if (d<5) m = 3000 else m = 10000
runif_cover = rvmf(m, rep(1,d), 0)
# Kovarianzmatrix
if (d<5) m = 1000 else m = 5000
runif_cover = runif_sphere(m, d)
# covariance matrix
Sigma_beta = rho_beta(runif_cover, runif_cover, beta)
Z_beta = mvrnorm(1, rep(0, m), Sigma_beta)
return(list("Realisierung" = max(Z_beta^2)))
return(list("realization" = max(Z_beta^2)))
}
#################################################################################################
## Berechnung kritischer Werte der Grenzverteilung von T per Monte Carlo ##
## calculation of critical values of the asymptotic distribution of T via Monte Carlo ##
# Umlaeufe
l = 10000
# Parameter
numCores = 2
# replication number l (cf. main manuscript, section 6) for critical values
l = 100000
# hyperparameters
numCores = 18
dim_grid = c(2, 3, 5, 10)
beta_grid = seq(1, 6)
grp_number = length(dim_grid)*length(beta_grid)
# Parameterkonstellationen fuer Monte Carlo
# parameter constellation for MonteCarlo package
parameter = list("d" = dim_grid, "beta" = beta_grid)
# Export benötigter Daten, Pakete und Funktionen
# export of required data, packages and functions
export_list = list("functions" = "rho_beta", "packages" = c("MASS", "Directional"))
# Monte Carlo Loops
# Monte Carlo loops
res_T_asy <- MonteCarlo::MonteCarlo(T_asymptotic_stat_value, nrep = l, param_list = parameter, ncpus = numCores, export_also = export_list)
# Erstellung eines Dataframes der Realisierungen
# producing a data.frame with all realizations
DF_T_asy <- MonteCarlo::MakeFrame(res_T_asy)
summary(DF_T_asy)
# In data.table umwandeln
# converting to a data.table
DT_T_asy <- setDT(DF_T_asy)
# Gruppen bilden
# creating groups
DT_T_asy[, grp := .GRP, by = c("d", "beta")]
DT_T_asy = DT_T_asy[order(grp, Realisierung)]
DT_T_asy = DT_T_asy[, list("Dimension" = d, "Beta" = beta, "Realisierung" = Realisierung,
"Quantil_90" = lapply(.SD, quantile, probs = 0.90),
"Quantil_95" = lapply(.SD, quantile, probs = 0.95),
"Quantil_99" = lapply(.SD, quantile, probs = 0.99)),
by = grp, .SDcols = c("Realisierung")]
DT_T_asy = DT_T_asy[order(grp, realization)]
DT_T_asy = DT_T_asy[, list("dimension" = d, "beta" = beta, "realization" = realization,
"level_10" = lapply(.SD, quantile, probs = 0.90),
"level_5" = lapply(.SD, quantile, probs = 0.95),
"level_1" = lapply(.SD, quantile, probs = 0.99)),
by = grp, .SDcols = c("realization")]
print(DT_T_asy)
save.image(file = "T_Statistic_GP_Daten.Rdata")
save.image(file = "T_Statistic_GP_data.Rdata")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment