-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
executable file
·121 lines (90 loc) · 3.61 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
---
output:
md_document:
variant: markdown_github
---
```{r, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE, comment = "#>", fig.path = "README/README-",
message = FALSE, warning = FALSE, fig.asp = 1, fig.align = 'center'
)
```
DirStats
========
```{r, echo = FALSE, results = 'asis'}
cat(
badger::badge_license(license = "GPLv3", color = "blue",
url = "https://www.gnu.org/licenses/gpl-3.0"),
badger::badge_github_actions(action = "R-CMD-check"),
badger::badge_cran_release(color = "green"),
badger::badge_cran_download(pkg = NULL, type = "grand-total"),
badger::badge_cran_download(pkg = NULL, type = "last-month")
)
```
<!-- <img src="" alt="DirStats hexlogo" align="right" width="200" style="padding: 0 15px; float: right;"/> -->
## Overview
Currently implementing nonparametric kernel density estimation, bandwidth selection, and other utilities for analyzing directional data. Further nonparametric tools expected to be included in subsequent releases.
## Installation
Get the released version from CRAN:
```{r, install-CRAN, eval = FALSE}
# Install the package
install.packages("DirStats")
# Load package
library(DirStats)
```
Alternatively, get the latest version from GitHub:
```{r, install-GitHub, eval = FALSE}
# Install the package
library(devtools)
install_github("egarpor/DirStats")
# Load package
library(DirStats)
```
```{r, load, echo = FALSE}
# Load package
library(DirStats)
```
## Usage
The following are examples of the usage of the Bai et al. (1988)'s kernel density estimator, the cross-validatory bandwidth selectors in Hall et al. (1987), and the plug-in bandwidth selectors in García-Portugués (2013).
### Compute bandwidths
```{r, bws}
# Sample from a von Mises--Fisher on S^2
q <- 2
n <- 300
set.seed(42)
samp <- rbind(rotasym::r_vMF(n = n / 3, mu = c(rep(0, q), 1), kappa = 5),
rotasym::r_vMF(n = n / 3, mu = c(rep(0, q), -1), kappa = 5),
rotasym::r_vMF(n = n / 3, mu = c(1, rep(0, q)), kappa = 5))
# LCV bandwidth
(h_LCV <- bw_dir_lcv(data = samp)$h_opt)
# LSCV bandwidth
(h_LSCV <- bw_dir_lscv(data = samp)$h_opt)
# ROT bandwidth
(h_ROT <- bw_dir_rot(data = samp))
# Mixture fit, for AMI and EMI bandwidth selectors
fit_mix <- bic_vmf_mix(data = samp)
# AMI bandwidth
(h_AMI <- bw_dir_ami(data = samp, fit_mix = fit_mix))
# EMI bandwidth
(h_EMI <- bw_dir_emi(data = samp, fit_mix = fit_mix)$h_opt)
```
### Compute kernel density estimator
```{r, kde, fig.asp = 2/3}
# Compute kde
l <- 200
th <- seq(0, 2 * pi, l = l)
phi <- seq(0, pi, l = l / 2)
ang <- expand.grid(th = th, phi = phi)
x <- to_sph(th = ang$th, ph = ang$phi)
kde <- kde_dir(x = x, data = samp, h = h_EMI)
# Visualization
contour(x = th, y = phi, z = matrix(kde, nrow = l, ncol = l / 2),
col = viridisLite::viridis(15),
levels = round(seq(min(kde), max(kde), length.out = 15), 4),
lwd = 2, xlab = expression(theta), ylab = expression(phi))
points(to_rad(samp), pch = 16)
```
## References
Bai, Z. D., Rao, C. R., and Zhao, L. C. (1988). Kernel estimators of density function of directional data. *Journal of Multivariate Analysis*, 27(1):24--39. [doi:10.1016/0047-259X(88)90113-3](https://doi.org/10.1016/0047-259X(88)90113-3).
García-Portugués, E. (2013). Exact risk improvement of bandwidth selectors for kernel density estimation with directional data. *Electronic Journal of Statistics*, 7:1655--1685. [doi:10.1214/13-ejs821](https://doi.org/10.1214/13-ejs821).
Hall, P., Watson, G. S., and Cabrera, J. (1987). Kernel density estimation with spherical data. *Biometrika*, 74(4):751--762. [doi:10.1093/biomet/74.4.751](https://doi.org/10.1093/biomet/74.4.751).