@@ -5,8 +5,8 @@ GPLv3](https://img.shields.io/badge/license-GPLv3-blue.svg)](https://www.gnu.org
5
5
[ ![ R build
6
6
status] ( https://github.com/egarpor/rotasym/workflows/R-CMD-check/badge.svg )] ( https://github.com/egarpor/rotasym/actions )
7
7
[ ![ ] ( https://www.r-pkg.org/badges/version/rotasym?color=green )] ( https://cran.r-project.org/package=rotasym )
8
- [ ![ ] ( http://cranlogs.r-pkg.org/badges/grand-total/rotasym?color=green )] ( https://cran.r-project.org/package=rotasym )
9
- [ ![ ] ( http://cranlogs.r-pkg.org/badges/last-month/rotasym?color=green )] ( https://cran.r-project.org/package=rotasym )
8
+ [ ![ ] ( http://cranlogs.r-pkg.org/badges/grand-total/rotasym )] ( https://cran.r-project.org/package=rotasym )
9
+ [ ![ ] ( http://cranlogs.r-pkg.org/badges/last-month/rotasym )] ( https://cran.r-project.org/package=rotasym )
10
10
11
11
<!-- <img src="" alt="rotasym hexlogo" align="right" width="200" style="padding: 0 15px; float: right;"/> -->
12
12
@@ -56,7 +56,7 @@ set.seed(123456789)
56
56
data_0 <- r_vMF(n = n , mu = theta , kappa = 1 )
57
57
```
58
58
59
- ### Specified-** θ ** case
59
+ ### Specified-θ case
60
60
61
61
``` r
62
62
# theta known
@@ -80,7 +80,7 @@ test_rotasym(data = data_0, theta = theta, type = "hyb_vMF")
80
80
# > Q_hyb_vMF = 46.329, df = 53, p-value = 0.7297
81
81
```
82
82
83
- ### Unspecified-** θ ** case
83
+ ### Unspecified-θ case
84
84
85
85
``` r
86
86
# theta unknown (employs the spherical mean as estimator)
@@ -107,7 +107,7 @@ test_rotasym(data = data_0, type = "hyb_vMF")
107
107
## Data application: test for the rotational symmetry of sunspots
108
108
109
109
The data application in García-Portugués, Paindaveine and Verdebout
110
- (2020) can be reproduced through the script
110
+ (2020) can be replicated through the script
111
111
[ sunspots-births.R] ( https://github.com/egarpor/rotasym/blob/master/data-raw/sunspots-births.R )
112
112
(data gathering and preprocessing) and the code snippet below.
113
113
@@ -183,12 +183,6 @@ example("sunspots_births")
183
183
# > snspt_+ radius = 1, type = "s", col = "lightblue", alpha = 0.25,
184
184
# > snspt_+ lit = FALSE)
185
185
# > snspt_+ }
186
- # > Error in dyn.load(dynlib <- getDynlib(dir)) :
187
- # > unable to load shared object '/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so':
188
- # > dlopen(/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so, 6): Symbol not found: _hb_buffer_add_utf8
189
- # > Referenced from: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so
190
- # > Expected in: flat namespace
191
- # > in /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so
192
186
# >
193
187
# > snspt_> n_cols <- 100
194
188
# >
@@ -199,12 +193,6 @@ example("sunspots_births")
199
193
# > snspt_> if (requireNamespace("rgl")) {
200
194
# > snspt_+ rgl::points3d(sunspots_23$X, col = viridisLite::viridis(n_cols)[cuts])
201
195
# > snspt_+ }
202
- # > Error in dyn.load(dynlib <- getDynlib(dir)) :
203
- # > unable to load shared object '/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so':
204
- # > dlopen(/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so, 6): Symbol not found: _hb_buffer_add_utf8
205
- # > Referenced from: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so
206
- # > Expected in: flat namespace
207
- # > in /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so
208
196
# >
209
197
# > snspt_> # Spörer's law: sunspots at the beginning of the solar cycle (dark blue
210
198
# > snspt_> # color) tend to appear at higher latitutes, gradually decreasing to the
@@ -218,76 +206,75 @@ example("sunspots_births")
218
206
# > snspt_> plot(kde <- density(x = V, bw = h, n = 2^13, from = -1, to = 1), col = 1,
219
207
# > snspt_+ xlim = c(-1, 1), ylim = c(0, 3), axes = FALSE, main = "",
220
208
# > snspt_+ xlab = "Cosines (latitude angles)", lwd = 2)
209
+ # >
210
+ # > snspt_> at <- seq(-1, 1, by = 0.25)
211
+ # >
212
+ # > snspt_> axis(2); axis(1, at = at)
213
+ # >
214
+ # > snspt_> axis(1, at = at, line = 1, tick = FALSE,
215
+ # > snspt_+ labels = paste0("(", 90 - round(acos(at) / pi * 180, 1), "º)"))
216
+ # >
217
+ # > snspt_> rug(V)
218
+ # >
219
+ # > snspt_> legend("topright", legend = c("Full cycle", "Initial 25% cycle",
220
+ # > snspt_+ "Final 25% cycle"),
221
+ # > snspt_+ lwd = 2, col = c(1, viridisLite::viridis(12)[c(3, 8)]))
222
+ # >
223
+ # > snspt_> # Density for the observations within the initial 25% of the cycle
224
+ # > snspt_> part1 <- sunspots_23$date < quantile(sunspots_23$date, 0.25)
225
+ # >
226
+ # > snspt_> V1 <- cosines(X = sunspots_23$X[part1, ], theta = c(0, 0, 1))
227
+ # >
228
+ # > snspt_> h1 <- bw.SJ(x = V1, method = "dpi")
229
+ # >
230
+ # > snspt_> lines(kde1 <- density(x = V1, bw = h1, n = 2^13, from = -1, to = 1),
231
+ # > snspt_+ col = viridisLite::viridis(12)[3], lwd = 2)
232
+ # >
233
+ # > snspt_> # Density for the observations within the final 25% of the cycle
234
+ # > snspt_> part2 <- sunspots_23$date > quantile(sunspots_23$date, 0.75)
235
+ # >
236
+ # > snspt_> V2 <- cosines(X = sunspots_23$X[part2, ], theta = c(0, 0, 1))
237
+ # >
238
+ # > snspt_> h2 <- bw.SJ(x = V2, method = "dpi")
239
+ # >
240
+ # > snspt_> lines(kde2 <- density(x = V2, bw = h2, n = 2^13, from = -1, to = 1),
241
+ # > snspt_+ col = viridisLite::viridis(12)[8], lwd = 2)
242
+ # >
243
+ # > snspt_> # Computation the level set of a kernel density estimator that contains
244
+ # > snspt_> # at least 1 - alpha of the probability (kde stands for an object
245
+ # > snspt_> # containing the output of density(x = data))
246
+ # > snspt_> kde_level_set <- function(kde, data, alpha) {
247
+ # > snspt_+
248
+ # > snspt_+ # Estimate c from alpha
249
+ # > snspt_+ c <- quantile(approx(x = kde$x, y = kde$y, xout = data)$y, probs = alpha)
250
+ # > snspt_+
251
+ # > snspt_+ # Begin and end index for the potentially many intervals in the level sets
252
+ # > snspt_+ kde_larger_c <- kde$y >= c
253
+ # > snspt_+ run_length_kde <- rle(kde_larger_c)
254
+ # > snspt_+ begin <- which(diff(kde_larger_c) > 0) + 1
255
+ # > snspt_+ end <- begin + run_length_kde$lengths[run_length_kde$values] - 1
256
+ # > snspt_+
257
+ # > snspt_+ # Return the [a_i, b_i], i = 1, ..., K in the K rows
258
+ # > snspt_+ return(cbind(kde$x[begin], kde$x[end]))
259
+ # > snspt_+
260
+ # > snspt_+ }
261
+ # >
262
+ # > snspt_> # Level set containing the 90% of the probability, in latitude angles
263
+ # > snspt_> 90 - acos(kde_level_set(kde = kde, data = V, alpha = 0.10)) / pi * 180
264
+ # > [,1] [,2]
265
+ # > [1,] -29.448244 -2.455986
266
+ # > [2,] 2.582017 28.123329
267
+ # >
268
+ # > snspt_> # Modes (in cosines and latitude angles)
269
+ # > snspt_> modes <- c(kde$x[kde$x < 0][which.max(kde$y[kde$x < 0])],
270
+ # > snspt_+ kde$x[kde$x > 0][which.max(kde$y[kde$x > 0])])
271
+ # >
272
+ # > snspt_> 90 - acos(modes) / pi * 180
273
+ # > [1] -13.69322 16.49001
221
274
```
222
275
223
276
<img src =" README/README-sunspots_births-1.png " style =" display : block ; margin : auto ;" />
224
277
225
- #>
226
- #> snspt_> at <- seq(-1, 1, by = 0.25)
227
- #>
228
- #> snspt_> axis(2); axis(1, at = at)
229
- #>
230
- #> snspt_> axis(1, at = at, line = 1, tick = FALSE,
231
- #> snspt_+ labels = paste0("(", 90 - round(acos(at) / pi * 180, 1), "º)"))
232
- #>
233
- #> snspt_> rug(V)
234
- #>
235
- #> snspt_> legend("topright", legend = c("Full cycle", "Initial 25% cycle",
236
- #> snspt_+ "Final 25% cycle"),
237
- #> snspt_+ lwd = 2, col = c(1, viridisLite::viridis(12)[c(3, 8)]))
238
- #>
239
- #> snspt_> # Density for the observations within the initial 25% of the cycle
240
- #> snspt_> part1 <- sunspots_23$date < quantile(sunspots_23$date, 0.25)
241
- #>
242
- #> snspt_> V1 <- cosines(X = sunspots_23$X[part1, ], theta = c(0, 0, 1))
243
- #>
244
- #> snspt_> h1 <- bw.SJ(x = V1, method = "dpi")
245
- #>
246
- #> snspt_> lines(kde1 <- density(x = V1, bw = h1, n = 2^13, from = -1, to = 1),
247
- #> snspt_+ col = viridisLite::viridis(12)[3], lwd = 2)
248
- #>
249
- #> snspt_> # Density for the observations within the final 25% of the cycle
250
- #> snspt_> part2 <- sunspots_23$date > quantile(sunspots_23$date, 0.75)
251
- #>
252
- #> snspt_> V2 <- cosines(X = sunspots_23$X[part2, ], theta = c(0, 0, 1))
253
- #>
254
- #> snspt_> h2 <- bw.SJ(x = V2, method = "dpi")
255
- #>
256
- #> snspt_> lines(kde2 <- density(x = V2, bw = h2, n = 2^13, from = -1, to = 1),
257
- #> snspt_+ col = viridisLite::viridis(12)[8], lwd = 2)
258
- #>
259
- #> snspt_> # Computation the level set of a kernel density estimator that contains
260
- #> snspt_> # at least 1 - alpha of the probability (kde stands for an object
261
- #> snspt_> # containing the output of density(x = data))
262
- #> snspt_> kde_level_set <- function(kde, data, alpha) {
263
- #> snspt_+
264
- #> snspt_+ # Estimate c from alpha
265
- #> snspt_+ c <- quantile(approx(x = kde$x, y = kde$y, xout = data)$y, probs = alpha)
266
- #> snspt_+
267
- #> snspt_+ # Begin and end index for the potentially many intervals in the level sets
268
- #> snspt_+ kde_larger_c <- kde$y >= c
269
- #> snspt_+ run_length_kde <- rle(kde_larger_c)
270
- #> snspt_+ begin <- which(diff(kde_larger_c) > 0)
271
- #> snspt_+ end <- begin + run_length_kde$lengths[run_length_kde$values] - 1
272
- #> snspt_+
273
- #> snspt_+ # Return the [a_i, b_i], i = 1, ..., K in the K rows
274
- #> snspt_+ return(cbind(kde$x[begin], kde$x[end]))
275
- #> snspt_+
276
- #> snspt_+ }
277
- #>
278
- #> snspt_> # Level set containing the 90% of the probability, in latitude angles
279
- #> snspt_> 90 - acos(kde_level_set(kde = kde, data = V, alpha = 0.10)) / pi * 180
280
- #> [,1] [,2]
281
- #> [1,] -29.464311 -2.469989
282
- #> [2,] 2.568013 28.107467
283
- #>
284
- #> snspt_> # Modes (in cosines and latitude angles)
285
- #> snspt_> modes <- c(kde$x[kde$x < 0][which.max(kde$y[kde$x < 0])],
286
- #> snspt_+ kde$x[kde$x > 0][which.max(kde$y[kde$x > 0])])
287
- #>
288
- #> snspt_> 90 - acos(modes) / pi * 180
289
- #> [1] -13.69322 16.49001
290
-
291
278
## References
292
279
293
280
García-Portugués, E., Paindaveine, D., and Verdebout, T. (2020). On
0 commit comments