op <- par(no.readonly = TRUE)
set.seed(1)
### generate a bimodal Gaussian mixture sample with 100000 points
n1 <- rbinom(1, 100000, .5)
x <- c(rnorm(n1), rnorm(100000-n1)/4+2)
# --------- Example 1: Grid evaluation --------------#
### evaluate exact and binned approximation on a
### grid and plot along with true density. We can
### use both Silverman's heuristic and the "mlcv" estimate.
xs <- seq(-5, 3.5, length = 1000)
ftrue <- dnorm(xs)/2 + dnorm(xs, 2, 1/4)/2
fhat <- fk_density(x, from = -5, to = 3.5)
fhat_bin <- fk_density(x, nbin = 1000, from = -5, to = 3.5)
par(mfrow = c(2, 2))
plot(xs, ftrue, type = 'l', lwd = 4, col = rgb(.7, .7, .7),
xlab = 'x', ylab = 'f(x)',
main = 'Exact evaluation, Silverman bandwidth')
lines(fhat, lty = 2)
plot(xs, ftrue, type = 'l', lwd = 4, col = rgb(.7, .7, .7),
xlab = 'x', ylab = 'f(x)',
main = 'Binned approximation, Silverman bandwidth')
lines(fhat_bin, lty = 2)
fhat <- fk_density(x, from = -5, to = 3.5, h = 'mlcv')
fhat_bin <- fk_density(x, nbin = 1000, from = -5,
to = 3.5, h = 'mlcv')
plot(xs, ftrue, type = 'l', lwd = 4, col = rgb(.7, .7, .7),
xlab = 'x', ylab = 'f(x)',
main = 'Exact evaluation, MLCV bandwidth')
lines(fhat, lty = 2)
plot(xs, ftrue, type = 'l', lwd = 4, col = rgb(.7, .7, .7),
xlab = 'x', ylab = 'f(x)',
main = 'Binned approximation, MLCV bandwidth')
lines(fhat_bin, lty = 2)
par(op)
# --------- Example 2: Evaluation at sample --------------#
### evaluate exact and binned approximation at the sample.
### Note that the output will be in the order of the original
### sample, and also the number of points will be large. It is
### not advisable, therefore, to simply plot these. We instead
### compute the mean squared deviations from the true density
ftrue <- sapply(x, function(xi) dnorm(xi)/2 + dnorm(xi, 2, 1/4)/2)
fhat <- fk_density(x, ngrid = NULL)
fhat_bin <- fk_density(x, nbin = 1000, x_eval = x)
mean((ftrue-fhat$y)^2)
mean((ftrue-fhat_bin$y)^2)
### now for MLCV bandwidth
fhat <- fk_density(x, h = 'mlcv', ngrid = NULL)
fhat_bin <- fk_density(x, nbin = 1000,
h = 'mlcv', x_eval = x)
mean((ftrue-fhat$y)^2)
mean((ftrue-fhat_bin$y)^2)
Run the code above in your browser using DataLab