From 20faf8c23a9ef122d460efdc428bfbdc7e63ce78 Mon Sep 17 00:00:00 2001
From: Jaroslav Borodavka <df2994@kit.edu>
Date: Fri, 15 Mar 2024 13:03:17 +0100
Subject: [PATCH] Adjusted code after addition of Bingham test
---
R/Comp_Statistics_Powertests_d_equal_2.R | 66 ++++++++++++------------
1 file changed, 34 insertions(+), 32 deletions(-)
diff --git a/R/Comp_Statistics_Powertests_d_equal_2.R b/R/Comp_Statistics_Powertests_d_equal_2.R
index 8907b7a..4d7ee4a 100644
--- a/R/Comp_Statistics_Powertests_d_equal_2.R
+++ b/R/Comp_Statistics_Powertests_d_equal_2.R
@@ -112,26 +112,26 @@ Bingham <- function(data){
return(n*d*(d+2)/2*(sum(diag(S %*% S)) - 1/d))
}
-# Alle Teststatistiken fuer MonteCarlo
+# collecting all test statistic in a list for MonteCarlo package
ls <- list(Kuiper, Watson, Ajne, Rayleigh, Cuesta_Albertos_multi, Bingham)
CircComp_stat_value <- function(data, index){
return(as.numeric(ls[[index]](data)))
}
-# Teststatistik mit gleichverteilten Daten
+# test realizations for uniform distribution
CircComp_unif <- function(n, dim, index){
return(list("Realisierung" = CircComp_stat_value(data = runif_sphere(n, dim),
index = index)))
}
-# Alternativen gemaeß Paper von Cutting, Pandaveine und Verdebout
+# alternatives according to paper by Cutting, Pandaveine und Verdebout ("TESTING UNIFORMITY ON HIGH-DIMENSIONAL SPHERES AGAINST MONOTONE ROTATIONALLY SYMMETRIC ALTERNATIVES")
CircComp_vMF <- function(n, dim, kappa, index){
return(list("Realisierung" = CircComp_stat_value(data = rvmf(n, c(1,rep(0, dim-1)), kappa),
index = index)))
}
-# Andere Alternativen
+# other alternatives
CircComp_mixvmf_1 <- function(n, dim, p, index){
U = runif(n, 0, 1)
sample = (U<p)*rvmf(n, c(-1,rep(0,dim-1)), 1) + (U>=p)*rvmf(n, c(1,rep(0,dim-1)), 1)
@@ -173,34 +173,34 @@ CircComp_LP <- function(n, dim, m, kappa, index){
}
#################################################################################################
-## Berechnung kritischer Werte per Monte Carlo unter der Nullhypothese ##
+## Calculation of critical values under the null hypothesis ##
-# Umlaeufe fuer kritische Werte
+# replication number l (cf. main manuscript) for critical values
reps = 20000
-# Parameter
+# hyperparameters
numCores = 18
n_grid = c(20, 50, 100, 500)
dim_grid = c(2)
-index_grid = seq(1, 5)
+index_grid = seq(1, 6)
grp_number = length(n_grid)*length(dim_grid)*length(index_grid)
-# Parameterkonstellationen fuer Monte Carlo
+# parameter constellation for MonteCarlo package
parameter = list("n" = n_grid, "dim" = dim_grid, "index" = index_grid)
-# Export benötigter Daten, Pakete und Funktionen
+# export of required data, packages and functions
export_list = list("functions" = c("euclidtoangular", "F_0", "Cuesta_Albertos_unique", "runif_sphere"),
"packages" = "MASS",
"data" = "rdirections")
# Monte Carlo Loops
res_CircComp_unif <- MonteCarlo::MonteCarlo(CircComp_unif, nrep = reps, param_list = parameter, ncpus = numCores, export_also = export_list)
-# Erstellung eines Dataframes der Realisierungen
+# producing a data.frame with all realizations
DF_CircComp_unif <- MonteCarlo::MakeFrame(res_CircComp_unif)
summary(DF_CircComp_unif)
-# In data.table umwandeln
+# converting to a data.table
DT_CircComp_unif <- setDT(DF_CircComp_unif)
-# Gruppen bilden
+# creating groups
DT_CircComp_unif[, grp := .GRP, by = c("n", "dim", "index")]
DT_CircComp_unif = DT_CircComp_unif[order(grp, Realisierung)]
DT_CircComp_unif_CA = DT_CircComp_unif[index == 5, list(n, "Dimension" = dim, "Index" = index, "Realisierung" = Realisierung,
@@ -208,7 +208,7 @@ DT_CircComp_unif_CA = DT_CircComp_unif[index == 5, list(n, "Dimension" = dim, "I
"Niveau_5" = lapply(.SD, quantile, probs = 0.05, na.rm = T),
"Niveau_1" = lapply(.SD, quantile, probs = 0.01, na.rm = T)),
by = grp, .SDcols = c("Realisierung")]
-DT_CircComp_unif = DT_CircComp_unif[index %in% seq(1,4), list(n, "Dimension" = dim, "Index" = index, "Realisierung" = Realisierung,
+DT_CircComp_unif = DT_CircComp_unif[index %in% c(seq(1,4), 6), list(n, "Dimension" = dim, "Index" = index, "Realisierung" = Realisierung,
"Niveau_10" = lapply(.SD, quantile, probs = 0.90, na.rm = T),
"Niveau_5" = lapply(.SD, quantile, probs = 0.95, na.rm = T),
"Niveau_1" = lapply(.SD, quantile, probs = 0.99, na.rm = T)),
@@ -218,9 +218,11 @@ print(DT_CircComp_unif)
#################################################################################################
-## Niveau-Test fuer konkurrierende Teststatistiken fuer d = 2 ##
+## level test for competing tests for d = 2 ##
-# Umlaeufe fuer Powertests
+# uniform distribution
+
+# replication number for power tests
a = 5000
CircComp_Test <- MonteCarlo::MonteCarlo(CircComp_unif, nrep = a, param_list = parameter, ncpus = numCores, export_also = export_list)
@@ -258,7 +260,7 @@ for (k in 1:grp_number) {
Power_CircComp[k, 2] = n
Power_CircComp[k, 3] = dim
- Power_CircComp[k, 4] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos")
+ Power_CircComp[k, 4] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos", "Bingham")
Power_CircComp[k, 5] = sum(dec_10[dec_10==1])/length(dec_10)
Power_CircComp[k, 6] = sum(dec_5[dec_5==1])/length(dec_5)
Power_CircComp[k, 7] = sum(dec_1[dec_1==1])/length(dec_1)
@@ -269,9 +271,9 @@ Power_CircComp = data.frame("Gruppe" = Power_CircComp[,1], "n" = Power_CircComp[
print(Power_CircComp)
#################################################################################################
-## Power-Test fuer konkurrierende Teststatistiken fuer d = 2 ##
+## power tests for competing tests for d = 2 ##
-# vMF
+# vMF alternative
kappa_grid = c(0.05, 0.1, 0.25, 0.5, 0.75, 1, 2)
grp_number = length(n_grid)*length(dim_grid)*length(kappa_grid)*length(index_grid)
@@ -316,7 +318,7 @@ for (k in 1:grp_number) {
Power_CircComp_vMF[k, 2] = umfang
Power_CircComp_vMF[k, 3] = dim
Power_CircComp_vMF[k, 4] = kappa
- Power_CircComp_vMF[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos")
+ Power_CircComp_vMF[k, 4] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos", "Bingham")
Power_CircComp_vMF[k, 6] = sum(dec_10[dec_10==1])/length(dec_10)
Power_CircComp_vMF[k, 7] = sum(dec_5[dec_5==1])/length(dec_5)
Power_CircComp_vMF[k, 8] = sum(dec_1[dec_1==1])/length(dec_1)
@@ -327,7 +329,7 @@ Power_CircComp_vMF = data.frame("Gruppe" = Power_CircComp_vMF[,1], "n" = Power_C
"Niveau_5" = Power_CircComp_vMF[,7], "Niveau_1" = Power_CircComp_vMF[,8])
print(Power_CircComp_vMF)
-# Mix FvML mit zwei Zentren (1)
+# Mix vMF with two centers (1)
p_grid = c(0.25, 0.5)
grp_number = length(n_grid)*length(dim_grid)*length(p_grid)*length(index_grid)
@@ -372,7 +374,7 @@ for (k in 1:grp_number) {
Power_CircComp_mixvmf_1[k, 2] = umfang
Power_CircComp_mixvmf_1[k, 3] = dim
Power_CircComp_mixvmf_1[k, 4] = prob
- Power_CircComp_mixvmf_1[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos")
+ Power_CircComp_mixvmf_1[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos", "Bingham")
Power_CircComp_mixvmf_1[k, 6] = sum(dec_10[dec_10==1])/length(dec_10)
Power_CircComp_mixvmf_1[k, 7] = sum(dec_5[dec_5==1])/length(dec_5)
Power_CircComp_mixvmf_1[k, 8] = sum(dec_1[dec_1==1])/length(dec_1)
@@ -383,7 +385,7 @@ Power_CircComp_mixvmf_1 = data.frame("Gruppe" = Power_CircComp_mixvmf_1[,1], "n"
"Niveau_5" = Power_CircComp_mixvmf_1[,7], "Niveau_1" = Power_CircComp_mixvmf_1[,8])
print(Power_CircComp_mixvmf_1)
-# Mix FvML mit zwei Zentren (2)
+# Mix vMF with two centers (2)
p_grid = c(0.5, 0.75)
grp_number = length(n_grid)*length(dim_grid)*length(p_grid)*length(index_grid)
@@ -428,7 +430,7 @@ for (k in 1:grp_number) {
Power_CircComp_mixvmf_2[k, 2] = umfang
Power_CircComp_mixvmf_2[k, 3] = dim
Power_CircComp_mixvmf_2[k, 4] = prob
- Power_CircComp_mixvmf_2[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos")
+ Power_CircComp_mixvmf_2[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos", "Bingham")
Power_CircComp_mixvmf_2[k, 6] = sum(dec_10[dec_10==1])/length(dec_10)
Power_CircComp_mixvmf_2[k, 7] = sum(dec_5[dec_5==1])/length(dec_5)
Power_CircComp_mixvmf_2[k, 8] = sum(dec_1[dec_1==1])/length(dec_1)
@@ -439,7 +441,7 @@ Power_CircComp_mixvmf_2 = data.frame("Gruppe" = Power_CircComp_mixvmf_2[,1], "n"
"Niveau_5" = Power_CircComp_mixvmf_2[,7], "Niveau_1" = Power_CircComp_mixvmf_2[,8])
print(Power_CircComp_mixvmf_2)
-# Mix FvML mit drei Zentren (1)
+# Mix vMF with three centers (1)
grp_number = length(n_grid)*length(dim_grid)*length(index_grid)
@@ -482,7 +484,7 @@ for (k in 1:grp_number) {
Power_CircComp_mixvmf_3[k, 2] = umfang
Power_CircComp_mixvmf_3[k, 3] = dim
Power_CircComp_mixvmf_3[k, 4] = 0.25
- Power_CircComp_mixvmf_3[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos")
+ Power_CircComp_mixvmf_3[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos", "Bingham")
Power_CircComp_mixvmf_3[k, 6] = sum(dec_10[dec_10==1])/length(dec_10)
Power_CircComp_mixvmf_3[k, 7] = sum(dec_5[dec_5==1])/length(dec_5)
Power_CircComp_mixvmf_3[k, 8] = sum(dec_1[dec_1==1])/length(dec_1)
@@ -493,7 +495,7 @@ Power_CircComp_mixvmf_3 = data.frame("Gruppe" = Power_CircComp_mixvmf_3[,1], "n"
"Niveau_1" = Power_CircComp_mixvmf_3[,8])
print(Power_CircComp_mixvmf_3)
-# Mix FvML mit drei Zentren (2)
+# Mix vMF with three centers (2)
grp_number = length(n_grid)*length(dim_grid)*length(index_grid)
@@ -534,7 +536,7 @@ for (k in 1:grp_number) {
Power_CircComp_mixvmf_4[k, 2] = umfang
Power_CircComp_mixvmf_4[k, 3] = dim
Power_CircComp_mixvmf_4[k, 4] = round(1/3, digits = 2)
- Power_CircComp_mixvmf_4[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos")
+ Power_CircComp_mixvmf_4[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos", , "Bingham")
Power_CircComp_mixvmf_4[k, 6] = sum(dec_10[dec_10==1])/length(dec_10)
Power_CircComp_mixvmf_4[k, 7] = sum(dec_5[dec_5==1])/length(dec_5)
Power_CircComp_mixvmf_4[k, 8] = sum(dec_1[dec_1==1])/length(dec_1)
@@ -591,7 +593,7 @@ for (k in 1:grp_number) {
Power_CircComp_bing_1[k, 2] = umfang
Power_CircComp_bing_1[k, 3] = dim
Power_CircComp_bing_1[k, 4] = tau
- Power_CircComp_bing_1[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos")
+ Power_CircComp_bing_1[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos", "Bingham")
Power_CircComp_bing_1[k, 6] = sum(dec_10[dec_10==1])/length(dec_10)
Power_CircComp_bing_1[k, 7] = sum(dec_5[dec_5==1])/length(dec_5)
Power_CircComp_bing_1[k, 8] = sum(dec_1[dec_1==1])/length(dec_1)
@@ -647,7 +649,7 @@ for (k in 1:grp_number) {
Power_CircComp_bing_2[k, 2] = umfang
Power_CircComp_bing_2[k, 3] = dim
Power_CircComp_bing_2[k, 4] = tau
- Power_CircComp_bing_2[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos")
+ Power_CircComp_bing_2[k, 5] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos", "Bingham")
Power_CircComp_bing_2[k, 6] = sum(dec_10[dec_10==1])/length(dec_10)
Power_CircComp_bing_2[k, 7] = sum(dec_5[dec_5==1])/length(dec_5)
Power_CircComp_bing_2[k, 8] = sum(dec_1[dec_1==1])/length(dec_1)
@@ -658,7 +660,7 @@ Power_CircComp_bing_2 = data.frame("Gruppe" = Power_CircComp_bing_2[,1], "n" = P
"Niveau_5" = Power_CircComp_bing_2[,7], "Niveau_1" = Power_CircComp_bing_2[,8])
print(Power_CircComp_bing_2)
-# Legendre-Polynom
+# Legendre polynomial
m_grid = c(3, 4)
kappa_grid = c(0.1, 0.5, 1)
@@ -706,7 +708,7 @@ for (k in 1:grp_number) {
Power_CircComp_LP[k, 3] = dim
Power_CircComp_LP[k, 4] = m
Power_CircComp_LP[k, 5] = kappa
- Power_CircComp_LP[k, 6] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos")
+ Power_CircComp_LP[k, 6] = switch(index, "Kuiper", "Watson", "Ajne", "Rayleigh", "Cuesta-Albertos", "Bingham")
Power_CircComp_LP[k, 7] = sum(dec_10[dec_10==1])/length(dec_10)
Power_CircComp_LP[k, 8] = sum(dec_5[dec_5==1])/length(dec_5)
Power_CircComp_LP[k, 9] = sum(dec_1[dec_1==1])/length(dec_1)
--
GitLab