智慧物流>> rotasym>> 返回
项目作者： egarpor

项目描述 ：
Tests for rotational symmetry on the hypersphere. Software companion for "On optimal tests for rotational symmetry against new classes of hyperspherical distributions"
高级语言： R
项目地址: git://github.com/egarpor/rotasym.git
创建时间： 2019-05-10T08:05:40Z
项目社区：https://github.com/egarpor/rotasym
开源协议：GNU General Public License v3.0
下载


rotasym

License: GPLv3
R build status


Overview

Software companion for the paper “On optimal tests for rotational
symmetry against new classes of hyperspherical distributions”
(García-Portugués, Paindaveine and Verdebout, 2020). It implements the
proposed tests for rotational symmetry of hyperspherical data and allows
to replicate the data application presented.

Installation

Get the released version from CRAN:

  1. # Install the package
  2. install.packages("rotasym")
  4. # Load package
  5. library(rotasym)

Alternatively, get the latest version from GitHub:

  1. # Install the package
  2. library(devtools)
  3. install_github("egarpor/rotasym")
  5. # Load package
  6. library(rotasym)

Tests usage

The following are some simple examples of the usage of the main function
of the package, test_rotasym, with simulated data. More examples are
available in ?test_rotasym.

  1. # Sample data from a vMF (rotational symmetric distribution about mu)
  2. n <- 200
  3. p <- 10
  4. theta <- c(1, rep(0, p - 1))
  5. set.seed(123456789)
  6. data_0 <- r_vMF(n = n, mu = theta, kappa = 1)

Specified-θ case

  1. # theta known
  2. test_rotasym(data = data_0, theta = theta, type = "sc")
  3. #> 
  4. #>  Scatter test for rotational symmetry
  5. #> 
  6. #> data:  data_0
  7. #> Q_sc = 35.013, df = 44, p-value = 0.8315
  8. test_rotasym(data = data_0, theta = theta, type = "loc_vMF")
  9. #> 
  10. #>  Location vMF test for rotational symmetry
  11. #> 
  12. #> data:  data_0
  13. #> Q_loc_vMF = 11.316, df = 9, p-value = 0.2547
  14. test_rotasym(data = data_0, theta = theta, type = "hyb_vMF")
  15. #> 
  16. #>  Hybrid vMF test (addition of statistics) for rotational symmetry
  17. #> 
  18. #> data:  data_0
  19. #> Q_hyb_vMF = 46.329, df = 53, p-value = 0.7297

Unspecified-θ case

  1. # theta unknown (employs the spherical mean as estimator)
  2. test_rotasym(data = data_0, type = "sc")
  3. #> 
  4. #>  Scatter test for rotational symmetry
  5. #> 
  6. #> data:  data_0
  7. #> Q_sc = 36.568, df = 44, p-value = 0.7793
  8. test_rotasym(data = data_0, type = "loc_vMF")
  9. #> 
  10. #>  Location vMF test for rotational symmetry
  11. #> 
  12. #> data:  data_0
  13. #> Q_loc_vMF = 12.335, df = 9, p-value = 0.1951
  14. test_rotasym(data = data_0, type = "hyb_vMF")
  15. #> 
  16. #>  Hybrid vMF test (addition of statistics) for rotational symmetry
  17. #> 
  18. #> data:  data_0
  19. #> Q_hyb_vMF = 48.902, df = 53, p-value = 0.6344

Data application: test for the rotational symmetry of sunspots

The data application in García-Portugués, Paindaveine and Verdebout
(2020) can be replicated through the script
sunspots-births.R
(data gathering and preprocessing) and the code snippet below.

  1. # Load data
  2. data("sunspots_births")
  3. sunspots_births$X <-
  4.   cbind(cos(sunspots_births$phi) * cos(sunspots_births$theta),
  5.         cos(sunspots_births$phi) * sin(sunspots_births$theta),
  6.         sin(sunspots_births$phi))
  8. # Test rotational symmetry for the 23rd cycle
  9. sunspots_23 <- subset(sunspots_births, cycle == 23)
  10. test_rotasym(data = sunspots_23$X, type = "sc", theta = c(0, 0, 1))
  11. #> 
  12. #>  Scatter test for rotational symmetry
  13. #> 
  14. #> data:  sunspots_23$X
  15. #> Q_sc = 3.5964, df = 2, p-value = 0.1656
  16. test_rotasym(data = sunspots_23$X, type = "loc", theta = c(0, 0, 1))
  17. #> 
  18. #>  Location test for rotational symmetry
  19. #> 
  20. #> data:  sunspots_23$X
  21. #> Q_loc = 1.5657, df = 2, p-value = 0.4571
  22. test_rotasym(data = sunspots_23$X, type = "hyb", theta = c(0, 0, 1))
  23. #> 
  24. #>  Hybrid test (addition of statistics) for rotational symmetry
  25. #> 
  26. #> data:  sunspots_23$X
  27. #> Q_hyb = 5.1622, df = 4, p-value = 0.2711
  29. # Test rotational symmetry for the 22nd cycle
  30. sunspots_22 <- subset(sunspots_births, cycle == 22)
  31. test_rotasym(data = sunspots_22$X, type = "sc", theta = c(0, 0, 1))
  32. #> 
  33. #>  Scatter test for rotational symmetry
  34. #> 
  35. #> data:  sunspots_22$X
  36. #> Q_sc = 4.4577, df = 2, p-value = 0.1077
  37. test_rotasym(data = sunspots_22$X, type = "loc", theta = c(0, 0, 1))
  38. #> 
  39. #>  Location test for rotational symmetry
  40. #> 
  41. #> data:  sunspots_22$X
  42. #> Q_loc = 8.7579, df = 2, p-value = 0.01254
  43. test_rotasym(data = sunspots_22$X, type = "hyb", theta = c(0, 0, 1))
  44. #> 
  45. #>  Hybrid test (addition of statistics) for rotational symmetry
  46. #> 
  47. #> data:  sunspots_22$X
  48. #> Q_hyb = 13.216, df = 4, p-value = 0.01027
  50. # More analyses in ?sunspots_births
  51. example("sunspots_births")
  52. #> 
  53. #> snspt_> # Load data
  54. #> snspt_> data("sunspots_births")
  55. #> 
  56. #> snspt_> # Transform to Cartesian coordinates
  57. #> snspt_> sunspots_births$X <-
  58. #> snspt_+   cbind(cos(sunspots_births$phi) * cos(sunspots_births$theta),
  59. #> snspt_+         cos(sunspots_births$phi) * sin(sunspots_births$theta),
  60. #> snspt_+         sin(sunspots_births$phi))
  61. #> 
  62. #> snspt_> # Plot data associated to the 23rd cycle
  63. #> snspt_> sunspots_23 <- subset(sunspots_births, cycle == 23)
  64. #> 
  65. #> snspt_> n <- nrow(sunspots_23$X)
  66. #> 
  67. #> snspt_> if (requireNamespace("rgl")) {
  68. #> snspt_+   rgl::plot3d(0, 0, 0, xlim = c(-1, 1), ylim = c(-1, 1), zlim = c(-1, 1),
  69. #> snspt_+               radius = 1, type = "s", col = "lightblue", alpha = 0.25,
  70. #> snspt_+               lit = FALSE)
  71. #> snspt_+ }
  72. #> 
  73. #> snspt_> n_cols <- 100
  74. #> 
  75. #> snspt_> cuts <- cut(x = sunspots_23$date, include.lowest = TRUE,
  76. #> snspt_+             breaks = quantile(sunspots_23$date,
  77. #> snspt_+                              probs = seq(0, 1, l = n_cols + 1)))
  78. #> 
  79. #> snspt_> if (requireNamespace("rgl")) {
  80. #> snspt_+   rgl::points3d(sunspots_23$X, col = viridisLite::viridis(n_cols)[cuts])
  81. #> snspt_+ }
  82. #> 
  83. #> snspt_> # Spörer's law: sunspots at the beginning of the solar cycle (dark blue
  84. #> snspt_> # color) tend to appear at higher latitudes, gradually decreasing to the
  85. #> snspt_> # equator as the solar cycle advances (yellow color)
  86. #> snspt_> 
  87. #> snspt_> # Estimation of the density of the cosines
  88. #> snspt_> V <- cosines(X = sunspots_23$X, theta = c(0, 0, 1))
  89. #> 
  90. #> snspt_> h <- bw.SJ(x = V, method = "dpi")
  91. #> 
  92. #> snspt_> plot(kde <- density(x = V, bw = h, n = 2^13, from = -1, to = 1), col = 1,
  93. #> snspt_+      xlim = c(-1, 1), ylim = c(0, 3), axes = FALSE, main = "",
  94. #> snspt_+      xlab = "Cosines (latitude angles)", lwd = 2)
  95. #> 
  96. #> snspt_> at <- seq(-1, 1, by = 0.25)
  97. #> 
  98. #> snspt_> axis(2); axis(1, at = at)
  99. #> 
  100. #> snspt_> axis(1, at = at, line = 1, tick = FALSE,
  101. #> snspt_+      labels = paste0("(", 90 - round(acos(at) / pi * 180, 1), "º)"))
  102. #> 
  103. #> snspt_> rug(V)
  104. #> 
  105. #> snspt_> legend("topright", legend = c("Full cycle", "Initial 25% cycle",
  106. #> snspt_+                               "Final 25% cycle"),
  107. #> snspt_+        lwd = 2, col = c(1, viridisLite::viridis(12)[c(3, 8)]))
  108. #> 
  109. #> snspt_> # Density for the observations within the initial 25% of the cycle
  110. #> snspt_> part1 <- sunspots_23$date < quantile(sunspots_23$date, 0.25)
  111. #> 
  112. #> snspt_> V1 <- cosines(X = sunspots_23$X[part1, ], theta = c(0, 0, 1))
  113. #> 
  114. #> snspt_> h1 <- bw.SJ(x = V1, method = "dpi")
  115. #> 
  116. #> snspt_> lines(kde1 <- density(x = V1, bw = h1, n = 2^13, from = -1, to = 1),
  117. #> snspt_+       col = viridisLite::viridis(12)[3], lwd = 2)
  118. #> 
  119. #> snspt_> # Density for the observations within the final 25% of the cycle
  120. #> snspt_> part2 <- sunspots_23$date > quantile(sunspots_23$date, 0.75)
  121. #> 
  122. #> snspt_> V2 <- cosines(X = sunspots_23$X[part2, ], theta = c(0, 0, 1))
  123. #> 
  124. #> snspt_> h2 <- bw.SJ(x = V2, method = "dpi")
  125. #> 
  126. #> snspt_> lines(kde2 <- density(x = V2, bw = h2, n = 2^13, from = -1, to = 1),
  127. #> snspt_+       col = viridisLite::viridis(12)[8], lwd = 2)
  128. #> 
  129. #> snspt_> # Computation the level set of a kernel density estimator that contains
  130. #> snspt_> # at least 1 - alpha of the probability (kde stands for an object
  131. #> snspt_> # containing the output of density(x = data))
  132. #> snspt_> kde_level_set <- function(kde, data, alpha) {
  133. #> snspt_+ 
  134. #> snspt_+   # Estimate c from alpha
  135. #> snspt_+   c <- quantile(approx(x = kde$x, y = kde$y, xout = data)$y, probs = alpha)
  136. #> snspt_+ 
  137. #> snspt_+   # Begin and end index for the potentially many intervals in the level sets
  138. #> snspt_+   kde_larger_c <- kde$y >= c
  139. #> snspt_+   run_length_kde <- rle(kde_larger_c)
  140. #> snspt_+   begin <- which(diff(kde_larger_c) > 0) + 1
  141. #> snspt_+   end <- begin + run_length_kde$lengths[run_length_kde$values] - 1
  142. #> snspt_+ 
  143. #> snspt_+   # Return the [a_i, b_i], i = 1, ..., K in the K rows
  144. #> snspt_+   return(cbind(kde$x[begin], kde$x[end]))
  145. #> snspt_+ 
  146. #> snspt_+ }
  147. #> 
  148. #> snspt_> # Level set containing the 90% of the probability, in latitude angles
  149. #> snspt_> 90 - acos(kde_level_set(kde = kde, data = V, alpha = 0.10)) / pi * 180
  150. #>            [,1]      [,2]
  151. #> [1,] -29.448244 -2.455986
  152. #> [2,]   2.582017 28.123329
  153. #> 
  154. #> snspt_> # Modes (in cosines and latitude angles)
  155. #> snspt_> modes <- c(kde$x[kde$x < 0][which.max(kde$y[kde$x < 0])],
  156. #> snspt_+            kde$x[kde$x > 0][which.max(kde$y[kde$x > 0])])
  157. #> 
  158. #> snspt_> 90 - acos(modes) / pi * 180
  159. #> [1] -13.69322  16.49001

References

García-Portugués, E., Paindaveine, D., and Verdebout, T. (2020). On
optimal tests for rotational symmetry against new classes of
hyperspherical distributions. Journal of the American Statistical
Association, 115(532):1873–1887.
doi:10.1080/01621459.2019.1665527.