From 61dd6d399495486e035aebaddc7b55fef5eef46f Mon Sep 17 00:00:00 2001
From: Jaroslav Borodavka <df2994@kit.edu>
Date: Wed, 17 Apr 2024 09:47:01 +0200
Subject: [PATCH] Superfluous script, hence removed.
---
R/Powertest_LP_Alternative.R | 266 -----------------------------------
1 file changed, 266 deletions(-)
delete mode 100644 R/Powertest_LP_Alternative.R
diff --git a/R/Powertest_LP_Alternative.R b/R/Powertest_LP_Alternative.R
deleted file mode 100644
index 252f929..0000000
--- a/R/Powertest_LP_Alternative.R
+++ /dev/null
@@ -1,266 +0,0 @@
-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")
-
--
GitLab