1#' @param bw The smoothing bandwidth to be used.
2#'   If numeric, the standard deviation of the smoothing kernel.
3#'   If character, a rule to choose the bandwidth, as listed in
4#'   [stats::bw.nrd()].
5#' @param adjust A multiplicate bandwidth adjustment. This makes it possible
6#'    to adjust the bandwidth while still using the a bandwidth estimator.
7#'    For example, `adjust = 1/2` means use half of the default bandwidth.
8#' @param kernel Kernel. See list of available kernels in [density()].
9#' @param n number of equally spaced points at which the density is to be
10#'   estimated, should be a power of two, see [density()] for
11#'   details
12#' @param trim If `FALSE`, the default, each density is computed on the
13#'   full range of the data. If `TRUE`, each density is computed over the
14#'   range of that group: this typically means the estimated x values will
15#'   not line-up, and hence you won't be able to stack density values.
16#'   This parameter only matters if you are displaying multiple densities in
17#'   one plot or if you are manually adjusting the scale limits.
18#' @section Computed variables:
19#' \describe{
20#'   \item{density}{density estimate}
21#'   \item{count}{density * number of points - useful for stacked density
22#'      plots}
23#'   \item{scaled}{density estimate, scaled to maximum of 1}
24#'   \item{ndensity}{alias for `scaled`, to mirror the syntax of
25#'    [`stat_bin()`]}
26#' }
27#' @export
28#' @rdname geom_density
29stat_density <- function(mapping = NULL, data = NULL,
30                         geom = "area", position = "stack",
31                         ...,
32                         bw = "nrd0",
33                         adjust = 1,
34                         kernel = "gaussian",
35                         n = 512,
36                         trim = FALSE,
37                         na.rm = FALSE,
38                         orientation = NA,
39                         show.legend = NA,
40                         inherit.aes = TRUE) {
41
42  layer(
43    data = data,
44    mapping = mapping,
45    stat = StatDensity,
46    geom = geom,
47    position = position,
48    show.legend = show.legend,
49    inherit.aes = inherit.aes,
50    params = list(
51      bw = bw,
52      adjust = adjust,
53      kernel = kernel,
54      n = n,
55      trim = trim,
56      na.rm = na.rm,
57      orientation = orientation,
58      ...
59    )
60  )
61}
62
63#' @rdname ggplot2-ggproto
64#' @format NULL
65#' @usage NULL
66#' @export
67StatDensity <- ggproto("StatDensity", Stat,
68  required_aes = "x|y",
69
70  default_aes = aes(x = after_stat(density), y = after_stat(density), fill = NA, weight = NULL),
71
72  setup_params = function(data, params) {
73    params$flipped_aes <- has_flipped_aes(data, params, main_is_orthogonal = FALSE, main_is_continuous = TRUE)
74
75    has_x <- !(is.null(data$x) && is.null(params$x))
76    has_y <- !(is.null(data$y) && is.null(params$y))
77    if (!has_x && !has_y) {
78      abort("stat_density() requires an x or y aesthetic.")
79    }
80
81    params
82  },
83
84  extra_params = c("na.rm", "orientation"),
85
86  compute_group = function(data, scales, bw = "nrd0", adjust = 1, kernel = "gaussian",
87                           n = 512, trim = FALSE, na.rm = FALSE, flipped_aes = FALSE) {
88    data <- flip_data(data, flipped_aes)
89    if (trim) {
90      range <- range(data$x, na.rm = TRUE)
91    } else {
92      range <- scales[[flipped_names(flipped_aes)$x]]$dimension()
93    }
94
95    density <- compute_density(data$x, data$weight, from = range[1],
96      to = range[2], bw = bw, adjust = adjust, kernel = kernel, n = n)
97    density$flipped_aes <- flipped_aes
98    flip_data(density, flipped_aes)
99  }
100
101)
102
103compute_density <- function(x, w, from, to, bw = "nrd0", adjust = 1,
104                            kernel = "gaussian", n = 512) {
105  nx <- length(x)
106  if (is.null(w)) {
107    w <- rep(1 / nx, nx)
108  } else {
109    w <- w / sum(w)
110  }
111
112  # if less than 2 points return data frame of NAs and a warning
113  if (nx < 2) {
114    warn("Groups with fewer than two data points have been dropped.")
115    return(new_data_frame(list(
116      x = NA_real_,
117      density = NA_real_,
118      scaled = NA_real_,
119      ndensity = NA_real_,
120      count = NA_real_,
121      n = NA_integer_
122    ), n = 1))
123  }
124
125  dens <- stats::density(x, weights = w, bw = bw, adjust = adjust,
126    kernel = kernel, n = n, from = from, to = to)
127
128  new_data_frame(list(
129    x = dens$x,
130    density = dens$y,
131    scaled =  dens$y / max(dens$y, na.rm = TRUE),
132    ndensity = dens$y / max(dens$y, na.rm = TRUE),
133    count =   dens$y * nx,
134    n = nx
135  ), n = length(dens$x))
136}
137