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

Adjusted code after addition of Bingham test

parent 42c04fa0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
R/Comp_Statistics_Powertests_d_equal_2.R
+ 34
− 32
View file @ 20faf8c2
......@@ -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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Help | Imprint | Privacy policy | Accessibility | Contact