# NOT RUN {
# Load data
data("sunspots_births")
# Transform to Cartesian coordinates
sunspots_births$X <-
cbind(cos(sunspots_births$phi) * cos(sunspots_births$theta),
cos(sunspots_births$phi) * sin(sunspots_births$theta),
sin(sunspots_births$phi))
# Plot data associated to the 23rd cycle
sunspots_23 <- subset(sunspots_births, cycle == 23)
n <- nrow(sunspots_23$X)
rgl::plot3d(0, 0, 0, xlim = c(-1, 1), ylim = c(-1, 1), zlim = c(-1, 1),
radius = 1, type = "s", col = "lightblue", alpha = 0.25,
lit = FALSE)
n_cols <- 100
cuts <- cut(x = sunspots_23$date, include.lowest = TRUE,
breaks = quantile(sunspots_23$date,
probs = seq(0, 1, l = n_cols + 1)))
rgl::points3d(sunspots_23$X, col = viridisLite::viridis(n_cols)[cuts])
# Sp<U+00F6>rer's law: sunspots at the beginning of the solar cycle (dark blue
# color) tend to appear at higher latitutes, gradually decreasing to the
# equator as the solar cycle advances (yellow color)
# Estimation of the density of the cosines
V <- cosines(X = sunspots_23$X, theta = c(0, 0, 1))
h <- bw.SJ(x = V, method = "dpi")
plot(kde <- density(x = V, bw = h, n = 2^13, from = -1, to = 1), col = 1,
xlim = c(-1, 1), ylim = c(0, 3), axes = FALSE, main = "",
xlab = "Cosines (latitude angles)", lwd = 2)
at <- seq(-1, 1, by = 0.25)
axis(2); axis(1, at = at)
axis(1, at = at, line = 1, tick = FALSE,
labels = paste0("(", 90 - round(acos(at) / pi * 180, 1), "<U+00BA>)"))
rug(V)
legend("topright", legend = c("Full cycle", "Initial 25% cycle",
"Final 25% cycle"),
lwd = 2, col = c(1, viridisLite::viridis(12)[c(3, 8)]))
# Density for the observations within the initial 25% of the cycle
part1 <- sunspots_23$date < quantile(sunspots_23$date, 0.25)
V1 <- cosines(X = sunspots_23$X[part1, ], theta = c(0, 0, 1))
h1 <- bw.SJ(x = V1, method = "dpi")
lines(kde1 <- density(x = V1, bw = h1, n = 2^13, from = -1, to = 1),
col = viridisLite::viridis(12)[3], lwd = 2)
# Density for the observations within the final 25% of the cycle
part2 <- sunspots_23$date > quantile(sunspots_23$date, 0.75)
V2 <- cosines(X = sunspots_23$X[part2, ], theta = c(0, 0, 1))
h2 <- bw.SJ(x = V2, method = "dpi")
lines(kde2 <- density(x = V2, bw = h2, n = 2^13, from = -1, to = 1),
col = viridisLite::viridis(12)[8], lwd = 2)
# Computation the level set of a kernel density estimator that contains
# at least 1 - alpha of the probability (kde stands for an object
# containing the output of density(x = data))
kde_level_set <- function(kde, data, alpha) {
# Estimate c from alpha
c <- quantile(approx(x = kde$x, y = kde$y, xout = data)$y, probs = alpha)
# Begin and end index for the potentially many intervals in the level sets
kde_larger_c <- kde$y >= c
run_length_kde <- rle(kde_larger_c)
begin <- which(diff(kde_larger_c) > 0)
end <- begin + run_length_kde$lengths[run_length_kde$values] - 1
# Return the [a_i, b_i], i = 1, ..., K in the K rows
return(cbind(kde$x[begin], kde$x[end]))
}
# Level set containing the 90% of the probability, in latitude angles
90 - acos(kde_level_set(kde = kde, data = V, alpha = 0.10)) / pi * 180
# Modes (in cosines and latitude angles)
modes <- c(kde$x[kde$x < 0][which.max(kde$y[kde$x < 0])],
kde$x[kde$x > 0][which.max(kde$y[kde$x > 0])])
90 - acos(modes) / pi * 180
# }
Run the code above in your browser using DataCamp Workspace