From 42c04fa011d8f4eb58c9d627499f8073fe1c0a7c Mon Sep 17 00:00:00 2001
From: Jaroslav Borodavka <df2994@kit.edu>
Date: Fri, 15 Mar 2024 12:17:46 +0100
Subject: [PATCH] Added Bingham test.

---
 R/Comp_Statistics_Powertests_d_equal_2.R | 46 +++++++++++++++---------
 1 file changed, 29 insertions(+), 17 deletions(-)

diff --git a/R/Comp_Statistics_Powertests_d_equal_2.R b/R/Comp_Statistics_Powertests_d_equal_2.R
index 8403088..8907b7a 100644
--- a/R/Comp_Statistics_Powertests_d_equal_2.R
+++ b/R/Comp_Statistics_Powertests_d_equal_2.R
@@ -1,16 +1,17 @@
+library(EnvStats)
 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")
 
-source("/Users/faik/Desktop/Mathematik/Masterarbeit/Ansatz_2/R-Codes/Legendre_Polynomial_Alternative.R")
+source("/home/jaroslav-borodavka/Schreibtisch/Dissertation/Arbeit/Bernoulli Submission/codes/codes_git_project/R/Legendre_Polynomial_Alternative.R")
 
 #################################################################################################
-## Konkurrierende Teststatistiken ## 
+## competing test statistics ## 
 
-# Testentscheid 
+# test decision 
 Decision <- function(teststatistic_value, critical){              
   
   decision <- 0
@@ -23,7 +24,7 @@ Decision <- function(teststatistic_value, critical){
   return(decision)
 }
 
-# Testentscheid CA-Teststatistik
+# test decision CA test statistic
 Decision_CA <- function(teststatistic_value, critical){              
   
   decision <- 0
@@ -36,7 +37,7 @@ Decision_CA <- function(teststatistic_value, critical){
   return(decision)
 }
 
-# Gleichverteilung auf S^(d-1)
+# uniformly distributed vectors on S^(d-1)
 runif_sphere <- function(n, d){
   X = mvrnorm(n, rep(0,d), diag(1,d))
   if (n == 1){
@@ -53,22 +54,22 @@ runif_sphere <- function(n, d){
   return(U)
 }
 
-# Transformation von Einheitsvektoren in zirkulƤre Daten via BogenmaƟ (nur d = 2), siehe n-sphere auf Wiki
+# transformation of unit vectors into circular data via radian (only d = 2)
 euclidtoangular <- function(unit_data){
   return(matrix(data = atan2(unit_data[,2], unit_data[,1]), ncol = 1))
 }
 
-# Kuiper-Teststatistik basierend auf Directional-Paket (nur d = 2)
+# Kuiper test based on Directional package (only d = 2), data in form of unit vectors
 Kuiper <- function(data){
-  return(Directional::kuiper(euclidtoangular(data), rads = TRUE)[1])
+  return(Directional::kuiper(euclidtoangular(data), rads = TRUE)$statistic)
 }
 
-# Watson-Teststatistik basierend auf Directional-Paket (nur d = 2)
+# Watson test based on Directional package (only d = 2), data in form of unit vectors
 Watson <- function(data){
-  return(Directional::watson(euclidtoangular(data), rads = TRUE)[1])
+  return(Directional::watson(euclidtoangular(data), rads = TRUE)$statistic)
 }
 
-# Cuesta-Albertos-Teststatistik, nur d = 2
+# Cuesta-Albertos test, only d = 2
 F_0 <- function(t){
   return((pi - acos(t))/pi)
 }
@@ -78,7 +79,7 @@ Cuesta_Albertos_unique <- function(data, rdirection){
   return(ks.test(projected_data, "F_0")$p.value)
 }
 
-# 100 gleichverteilte Richtungen
+# 25 uniformly distributed directions
 rdirections = runif_sphere(25, 2)
 
 Cuesta_Albertos_multi <- function(data){
@@ -86,12 +87,12 @@ Cuesta_Albertos_multi <- function(data){
   return(min(p_values))
 }
 
-# Rayleigh-Teststatistik basierend auf Directional-Paket 
+# Rayleigh based on Directional package 
 Rayleigh <- function(data){
-  return(Directional::rayleigh(data, modif = TRUE)[1])
+  return(Directional::rayleigh(data, modif = TRUE)$statistic)
 }
 
-# Ajne-Teststatistik
+# Ajne test
 Ajne <- function(data){
   n = dim(data)[1]
   temp = data %*% t(data)
@@ -99,9 +100,20 @@ Ajne <- function(data){
   return(n/4 - 1/(n*pi)*sum(Psi))
 }
 
+# Bingham test 
+Bingham <- function(data){
+  n = dim(data)[1]
+  d = dim(data)[2]
+  temp = matrix(0, d, d)
+  for (i in 1:n) {
+    temp = temp + data[i,] %*% t(data[i,])
+  }
+  S = 1/n*temp
+  return(n*d*(d+2)/2*(sum(diag(S %*% S)) - 1/d))
+}
 
 # Alle Teststatistiken fuer MonteCarlo
-ls <-  list(Kuiper, Watson, Ajne, Rayleigh, Cuesta_Albertos_multi)
+ls <-  list(Kuiper, Watson, Ajne, Rayleigh, Cuesta_Albertos_multi, Bingham)
 
 CircComp_stat_value <- function(data, index){
   return(as.numeric(ls[[index]](data)))
-- 
GitLab