From fda26fb1ba9170d5abdf84a7f9a076ec0e244834 Mon Sep 17 00:00:00 2001
From: Jaroslav Borodavka <df2994@kit.edu>
Date: Sat, 13 Apr 2024 13:50:35 +0200
Subject: [PATCH] Translated to German and added a preamble.
---
R/T_Statistic_GP.R | 76 +++++++++++++++++++++++++++++++---------------
1 file changed, 51 insertions(+), 25 deletions(-)
diff --git a/R/T_Statistic_GP.R b/R/T_Statistic_GP.R
index 1985776..8867be2 100644
--- a/R/T_Statistic_GP.R
+++ b/R/T_Statistic_GP.R
@@ -1,14 +1,21 @@
+#################################################################################################
+#################################################################################################
+##
+#################################################################################################
+#################################################################################################
+
+# 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")
--
GitLab