Superfluous script, hence removed.

R/Powertest_LP_Alternative.R

-library(Directional)
-library(MASS)
-library(MonteCarlo)
-library(data.table)
-library(zipfR)        # incomplete Beta function Ibeta
-
-setwd("/Users/faik/Desktop/Mathematik/Masterarbeit/Ansatz_2/R-Codes")
-
-source("/Users/faik/Desktop/Mathematik/Masterarbeit/Ansatz_2/R-Codes/Legendre_Polynomial_Alternative.R")
-
-#################################################################################################
-## Funktionen und Teststatistik ## 
-
-# Erwartungswerte Psi_d
-psi_d <- function(d, beta){
-  return(if (beta%%2 == 0) gamma((beta+1)/2)*gamma(d/2)/(sqrt(pi)*gamma((beta+d)/2)) else 0) 
-}
-
-# Quadrierte und zentrierte Projektionen
-projection <- function(unif_point_cover, data, beta){
-  n = dim(data)[1]
-  d = dim(data)[2]
-  scalarproducts = data %*% t(unif_point_cover)
-  summen = colSums((scalarproducts)^beta)
-  return((1/n*summen - psi_d(d, beta))^2)
-}
-
-# Teststatistik T_n,beta von Ebner, cover_n ist die Anzahl gleichverteilter Punkte auf der Sphäre
-T_stat_value <- function(data, beta){
-  n = dim(data)[1]
-  d = dim(data)[2]
-  cover_n = 0
-  if (d < 5) cover_n = 5000 else cover_n = 20000
-  unif_points = rvmf(cover_n, rep(1,d), 0)
-  projections = projection(unif_point_cover = unif_points, data = data, beta = beta)
-  return(n*max(projections))
-}
-
-T_LP <- function(n, dim, m, kappa, beta){
-  return(list("Realisierung" = T_stat_value(data = rLP(n, m, c(1,rep(0,dim-1)), kappa/sqrt(n)), beta = beta)))
-}
-
-SpherComp_LP <- function(n, dim, m, kappa, index){
-  return(list("Realisierung" = SpherComp_stat_value(data = rLP(n, m, c(1,rep(0,dim-1)), kappa), index = index)))
-}
-
-# Ajne-Teststatistik
-Ajne <- function(data){
-  n = dim(data)[1]
-  temp = data %*% t(data)
-  Psi = acos(temp[upper.tri(temp)])
-  return(n/4 - 1/(n*pi)*sum(Psi))
-}
-
-# Rayleigh-Teststatistik basierend auf Directional-Paket 
-Rayleigh <- function(data){
-  return(Directional::rayleigh(data, modif = TRUE)[1])
-}
-
-# Bingham-Teststatistik
-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))
-}
-
-# Gine-Teststatistik G
-Gine_G <- function(data){
-  n = dim(data)[1]
-  d = dim(data)[2]
-  temp = data %*% t(data)
-  Psi = acos(temp[upper.tri(temp)])
-  return(n/2 - (d-1)*gamma((d-1)/2)^2/(2*n*gamma(d/2)^2)*sum(sin(Psi)))
-}
-
-# Cuesta-Albertos-Teststatistik
-F_CA <- function(d){
-  if (d==2){
-    res <- function(t) 1 - acos(t)/pi
-  } else if(d==3){
-    res <- function(t) (t+1)/2
-  } else{
-    res <- function(t) 1/2*(1 + sign(t)*Ibeta(t^2, a = 1/2, b = (d-1)/2)/Cbeta(a = 1/2, b = (d-1)/2))
-  }
-  return(res)
-}
-
-Cuesta_Albertos_unique <- function(data, rdirection){
-  d = dim(data)[2]
-  cdf = F_CA(d)
-  projected_data = data %*% rdirection
-  return(ks.test(projected_data, "cdf")$p.value)
-}
-
-# 100 gleichverteilte Richtungen fuer die CA-Teststatistik
-rdirections_2 <- runif_sphere(25, 2)
-rdirections_3 <- runif_sphere(100, 3)
-rdirections_5 <- runif_sphere(100, 5)
-rdirections_10 <- runif_sphere(100, 10)
-
-Cuesta_Albertos_multi <- function(data){
-  d = dim(data)[2]
-  if (d==2){
-    rdirections = rdirections_2
-  } else if (d==3){
-    rdirections = rdirections_3
-  } else if (d==5){
-    rdirections = rdirections_5
-  } else {
-    rdirections = rdirections_10
-  }
-  p_values = apply(rdirections, 1, Cuesta_Albertos_unique, data = data)
-  return(min(p_values))
-}
-
-# Cramér-von Mises-Teststatistik
-integrand <- function(t, angle, dim){
-  return(F_CA(dim)(t)*F_CA(dim-1)(t*tan(angle/2)/sqrt(1-t^2))/Cbeta(a = 1/2, b = (dim-1)/2)*(1-t^2)^((dim-3)/2))
-}
-
-psi_CvM <- function(d){
-  if (d==2){
-    res <- function(theta) 1/2 + theta/(2*pi)*(theta/(2*pi) - 1)
-  } else if(d==3){
-    res <- function(theta) 1/2 - 1/4*sin(theta/2)
-  } else{
-    res <- function(theta){
-      return(-3/4 + theta/(2*pi) + 2*(F_CA(d)(cos(theta/2)))^2 - 4*integrate(integrand, 
-                                                                             lower = 0, 
-                                                                             upper = cos(theta/2), 
-                                                                             angle = theta,
-                                                                             dim = d)$value)
-    }
-  }
-  return(res)
-}
-
-CramerVonMises <- function(data){
-  n = dim(data)[1]
-  d = dim(data)[2]
-  temp = data %*% t(data)
-  func = psi_CvM(d)
-  psi_q_acos = sapply(acos(temp[upper.tri(temp)]), func)
-  return(2/n*sum(psi_q_acos) + (3-2*n)/6)
-}
-
-# Transformation von Einheitsvektoren in zirkuläre Daten via Bogenmaß (nur d = 2), siehe n-sphere auf Wiki
-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 <- function(data){
-  return(Directional::kuiper(euclidtoangular(data), rads = TRUE)[1])
-}
-
-# Watson-Teststatistik basierend auf Directional-Paket (nur d = 2)
-Watson <- function(data){
-  return(Directional::watson(euclidtoangular(data), rads = TRUE)[1])
-}
-
-# Alle Teststatistiken fuer MonteCarlo
-ls_spher <-  list(Ajne, Rayleigh, Bingham, Gine_G, Cuesta_Albertos_multi, CramerVonMises)
-
-# Alle Teststatistiken fuer MonteCarlo
-ls_circ <-  list(Kuiper, Watson, Ajne, Rayleigh, Cuesta_Albertos_multi)
-
-SpherComp_stat_value <- function(data, index){
-  return(as.numeric(ls_spher[[index]](data)))
-}
-
-CircComp_stat_value <- function(data, index){
-  return(as.numeric(ls_circ[[index]](data)))
-}
-
-T_LP <- function(n, dim, m, kappa, beta){
-  return(list("Realisierung" = T_stat_value(data = rLP(n, m, c(1,rep(0,dim-1)), kappa), beta = beta)))
-}
-
-SpherComp_LP <- function(n, dim, m, kappa, index){
-  return(list("Realisierung" = SpherComp_stat_value(data = rLP(n, m, c(1,rep(0,dim-1)), kappa), index = index)))
-}
-
-CircComp_LP <- function(n, dim, m, kappa, index){
-  return(list("Realisierung" = CircComp_stat_value(data = rLP(n, m, c(1,rep(0,dim-1)), kappa), index = index)))
-}
-
-#################################################################################################
-## Berechnung kritischer Werte per Monte Carlo unter der Nullhypothese ## 
-
-# Parameter
-a = 5000
-numCores = 2
-n_grid = c(20, 50, 100, 500)                              
-dim_grid = c(2, 3, 5, 10)
-beta_grid = seq(1, 6)
-m_grid = seq(2, 6)
-kappa_grid = c(0.1, 0.5, 1)
-
-parameter = list("n" = n_grid, "dim" = dim_grid, "m" = m_grid, "kappa" = kappa_grid, "beta" = beta_grid)    
-
-export_list = list("functions" = c("psi_d", "projection", "euclidtoangular", "F_CA", "Cuesta_Albertos_unique", "integrand", "psi_CvM"),
-                   "packages" = c("MASS", "Directional",  "zipfR"),
-                   "data" = c("rdirections_2", "rdirections_3", "rdirections_5", "rdirections_10", "ls_circ", "ls_spher"))
-
-# T_n,beta
-
-parameter = list("n" = n_grid, "dim" = dim_grid, "m" = m_grid, "kappa" = kappa_grid, "beta" = beta_grid)    
-
-res_T_LP_neu <- MonteCarlo::MonteCarlo(T_LP, nrep = a, param_list = parameter, 
-                                   ncpus = numCores, max_grid = 1500, time_n_test = TRUE, export_also = export_list)     
-DF_T_LP_neu <- MonteCarlo::MakeFrame(res_T_LP_neu)                
-summary(DF_T_LP_neu)
-
-DT_T_LP_neu <- setDT(DF_T_LP_neu)
-print(DT_T_LP_neu)
-
-DT_T_LP_neu[, grp := .GRP, by = c("n", "dim", "m", "kappa", "beta")]                
-DT_T_LP_neu = DT_T_LP_neu[order(grp)]
-DT_T_LP_neu = na.omit(DT_T_LP_neu)
-
-# Konkurrierende Teststatistiken, d > 2
-
-dim_grid = c(3, 5, 10)
-index_grid = seq(1, 6)
-
-parameter = list("n" = n_grid, "dim" = dim_grid, "m" = m_grid, "kappa" = kappa_grid, "index" = index_grid) 
-
-res_SpherComp_LP_neu <- MonteCarlo::MonteCarlo(SpherComp_LP, nrep = a, param_list = parameter, 
-                                           ncpus = numCores, max_grid = 1500, export_also = export_list)     
-DF_SpherComp_LP_neu <- MonteCarlo::MakeFrame(res_SpherComp_LP_neu)                
-summary(DF_SpherComp_LP_neu)
-
-DT_SpherComp_LP_neu <- setDT(DF_SpherComp_LP_neu)
-print(DT_SpherComp_LP_neu)
-
-DT_SpherComp_LP_neu[, grp := .GRP, by = c("n", "dim", "m", "kappa", "index")]                
-DT_SpherComp_LP_neu = DT_SpherComp_LP_neu[order(grp)]
-DT_SpherComp_LP_neu = na.omit(DT_SpherComp_LP_neu)
-
-# Konkurrierende Teststatistiken, d = 2
-
-dim_grid = c(2)
-index_grid = seq(1, 5)
-
-parameter = list("n" = n_grid, "dim" = dim_grid, "m" = m_grid, "kappa" = kappa_grid, "index" = index_grid) 
-
-res_CircComp_LP_neu <- MonteCarlo::MonteCarlo(CircComp_LP, nrep = a, param_list = parameter,
-                                          ncpus = numCores, export_also = export_list)     
-DF_CircComp_LP_neu <- MonteCarlo::MakeFrame(res_CircComp_LP_neu)                
-summary(DF_CircComp_LP_neu)
-
-DT_CircComp_LP_neu <- setDT(DF_CircComp_LP_neu)
-print(DT_CircComp_LP_neu)
-
-DT_CircComp_LP_neu[, grp := .GRP, by = c("n", "dim", "m", "kappa", "index")]                
-DT_CircComp_LP_neu = DT_CircComp_LP_neu[order(grp)]
-DT_CircComp_LP_neu = na.omit(DT_CircComp_LP_neu)
-
-save.image(file = "Powertest_LP_Alternative_Daten.RData")
-