rmvg {kernelboot}

Random generation from multivariate Gaussian kernel density

Description

Random generation from multivariate Gaussian kernel density

Usage

rmvg(n, y, bw = bw.silv(y), weights = NULL, adjust = 1)

Arguments

n

number of observations. If length(n) > 1, the length is taken to be the number required.

y

numeric matrix or data.frame.

bw

numeric matrix with number of rows and columns equal to ncol(y); the smoothing bandwidth to be used. This is the covariance matrix of the smoothing kernel. If provided as a single value, the same bandwidth is used for each variable. If provided as a single value, or as a vector, variables are considered as uncorrelated.

weights

numeric vector of length equal to nrow(y); must be non-negative.

adjust

scalar; the bandwidth used is actually adjust*bw. This makes it easy to specify values like 'half the default' bandwidth.

Details

Multivariate kernel density estimator with multivariate Gaussian (normal) kernels K_H is defined as

\hat{f_H}(\mathbf{x}) = \sum_{i=1}^n w_i \, K_H \left( \mathbf{x}-\boldsymbol{y}_i \right)

where w is a vector of weights such that all w_i \ge 0 and \sum_i w_i = 1 (by default uniform 1/n weights are used), K_H is kernel K parametrized by bandwidth matrix H and \boldsymbol{y} is a matrix of data points used for estimating the kernel density.

Random generation from multivariate normal distribution is possible by taking

x = A' z + \mu

where z is a vector of m i.i.d. standard normal deviates, \mu is a vector of means and A is a m \times m matrix such that A'A=\Sigma (A is a Cholesky factor of \Sigma). In the case of multivariate Gaussian kernel density, \mu, is the i-th row of \boldsymbol{y}, where i is drawn randomly with replacement with probability proportional to w_i, and \Sigma is the bandwidth matrix H.

For functions estimating kernel densities please check KernSmooth, ks, or other packages reviewed by Deng and Wickham (2011).

References

Deng, H. and Wickham, H. (2011). Density estimation in R. http://vita.had.co.nz/papers/density-estimation.pdf

See Also

kernelboot

Examples

set.seed(1)

dat <- mtcars[, c(1,3)]
bw <- bw.silv(dat)
X <- rmvg(5000, dat, bw = bw)

if (requireNamespace("ks", quietly = TRUE)) {

   pal <- colorRampPalette(c("chartreuse4", "yellow", "orange", "brown"))
   col <- pal(10)[cut(ks::kde(dat, H = bw, eval.points = X)$estimate, breaks = 10)]

   plot(X, col = col, pch = 19, axes = FALSE,
        main = "Multivariate Gaussian Kernel")
   points(dat, pch = 2, col = "blue")
   axis(1); axis(2)

} else {

   plot(X, pch = 16, axes = FALSE, col = "#458B004D",
        main = "Multivariate Gaussian Kernel")
   points(dat, pch = 2, col = "red", lwd = 2)
   axis(1); axis(2)

}