Density Estimation

This function is a wrapper over different methods of density estimation. By default, it uses the base R density with by default uses a different smoothing bandwidth ("SJ") from the legacy default implemented the base R density function ("nrd0"). However, Deng \& Wickham suggest that method = "KernSmooth" is the fastest and the most accurate.

estimate_density(
  x,
  method = "kernel",
  precision = 2^10,
  extend = FALSE,
  extend_scale = 0.1,
  bw = "SJ",
  ...
)

# S3 method for data.frame
estimate_density(
  x,
  method = "kernel",
  precision = 2^10,
  extend = FALSE,
  extend_scale = 0.1,
  bw = "SJ",
  group_by = NULL,
  ...
)

Arguments

x	Vector representing a posterior distribution, or a data frame of such vectors. Can also be a Bayesian model (stanreg, brmsfit, MCMCglmm, mcmc or bcplm) or a BayesFactor model.
method	Density estimation method. Can be "kernel" (default), "logspline" or "KernSmooth".
precision	Number of points of density data. See the n parameter in density.
extend	Extend the range of the x axis by a factor of extend_scale.
extend_scale	Ratio of range by which to extend the x axis. A value of 0.1 means that the x axis will be extended by 1/10 of the range of the data.
bw	See the eponymous argument in density. Here, the default has been changed for "SJ", which is recommended.
...	Currently not used.
group_by	Optional character vector. If not NULL and x is a data frame, density estimation is performed for each group (subset) indicated by group_by.

Note

There is also a plot()-method implemented in the see-package.

References

Deng, H., & Wickham, H. (2011). Density estimation in R. Electronic publication.

Examples

library(bayestestR)

set.seed(1)
x <- rnorm(250, 1)

# Methods
density_kernel <- estimate_density(x, method = "kernel")
density_logspline <- estimate_density(x, method = "logspline")
density_KernSmooth <- estimate_density(x, method = "KernSmooth")
#> Loading required namespace: KernSmooth
density_mixture <- estimate_density(x, method = "mixture")
#> Loading required namespace: mclust

hist(x, prob = TRUE)
lines(density_kernel$x, density_kernel$y, col = "black", lwd = 2)
lines(density_logspline$x, density_logspline$y, col = "red", lwd = 2)
lines(density_KernSmooth$x, density_KernSmooth$y, col = "blue", lwd = 2)
lines(density_mixture$x, density_mixture$y, col = "green", lwd = 2)

# Extension
density_extended <- estimate_density(x, extend = TRUE)
density_default <- estimate_density(x, extend = FALSE)

hist(x, prob = TRUE)
lines(density_extended$x, density_extended$y, col = "red", lwd = 3)
lines(density_default$x, density_default$y, col = "black", lwd = 3)

df <- data.frame(replicate(4, rnorm(100)))
head(estimate_density(df))
#>   Parameter         x          y
#> 1        X1 -2.403096 0.03581143
#> 2        X1 -2.398158 0.03603563
#> 3        X1 -2.393219 0.03625900
#> 4        X1 -2.388280 0.03648072
#> 5        X1 -2.383342 0.03670092
#> 6        X1 -2.378403 0.03692027
if (FALSE) {
# rstanarm models
# -----------------------------------------------
library(rstanarm)
model <- stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0)
head(estimate_density(model))

library(emmeans)
head(estimate_density(emtrends(model, ~1, "wt")))

# brms models
# -----------------------------------------------
library(brms)
model <- brms::brm(mpg ~ wt + cyl, data = mtcars)
estimate_density(model)
}