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