Introduction to epkde

Overview

epkde implements the approximate Bayesian method for bandwidth selection in multivariate kernel density estimation (KDE) from Filippone & Sanguinetti (2011). The key idea is to place a prior on the kernel precision matrix \(\Lambda\) (the inverse bandwidth), and use the Expectation Propagation (EP) algorithm to compute an approximate posterior, exploiting a leave-one-out cross-validated likelihood.

Three kernel structures are supported:

Structure	Prior on \(\Lambda\)	Free parameters
Isotropic	Gamma on scalar \(\lambda\)	1
Diagonal	Independent Gammas on \(\lambda_k\)	\(d\)
Full	Wishart on matrix \(\Lambda\)	\(d(d+1)/2\)

The model evidence returned by each fit makes it possible to select the appropriate structure using Bayes factors, without a held-out test set.

Quick start: bivariate Gaussian mixture

# Generate data from a bivariate Gaussian mixture
n <- 150; d <- 2
x <- rbind(matrix(rnorm(n/2 * d, mean = c(-2,  0)), ncol = d),
           matrix(rnorm(n/2 * d, mean = c( 2,  0)), ncol = d))

Isotropic fit

fit_iso <- ep_kde_isotropic(x, prior_shape = 1, prior_rate = 1)
cat("Posterior mean precision (lambda):", fit_iso$post_shape / fit_iso$post_rate, "\n")
#> Posterior mean precision (lambda): 2.316491
cat("Log model evidence:", fit_iso$log_evidence, "\n")
#> Log model evidence: -589.318

Diagonal fit

fit_diag <- ep_kde_diagonal(x,
                             prior_shape = rep(1, d),
                             prior_rate  = rep(1, d))
cat("Posterior mean precision per dimension:",
    fit_diag$post_shape / fit_diag$post_rate, "\n")
#> Posterior mean precision per dimension: 2.28161 2.224204
cat("Log model evidence:", fit_diag$log_evidence, "\n")
#> Log model evidence: -590.5628

Full precision matrix fit

fit_full <- ep_kde_full(x, prior_nu = 0)
cat("Posterior nu:", fit_full$post_nu, "\n")
#> Posterior nu: 131.6236
cat("Posterior mean precision matrix:\n")
#> Posterior mean precision matrix:
print(fit_full$post_mean)
#>               [,1]          [,2]
#> [1,]  0.0921679865 -0.0002349248
#> [2,] -0.0002349248  0.0906878646
cat("Log model evidence:", fit_full$log_evidence, "\n")
#> Log model evidence: -799.2478

Model comparison via Bayes factors

ev <- c(isotropic = model_evidence(fit_iso),
        diagonal  = model_evidence(fit_diag),
        full      = model_evidence(fit_full))
cat("Log evidences:\n"); print(ev)
#> Log evidences:
#> isotropic  diagonal      full 
#> -589.3180 -590.5628 -799.2478
cat("\nLog Bayes factor (diagonal vs. isotropic):",
    ev["diagonal"] - ev["isotropic"], "\n")
#> 
#> Log Bayes factor (diagonal vs. isotropic): -1.244776

Evaluating the fitted density

# Grid for evaluation
gr    <- seq(-6, 6, length.out = 50)
grid  <- as.matrix(expand.grid(gr, gr))

# Use the diagonal fit: Lambda = diag(post_shape / post_rate)
Lambda_diag <- diag(fit_diag$post_shape / fit_diag$post_rate)
p_hat       <- kde_predict(grid, x, Lambda_diag)

# Contour plot
contour(gr, gr, matrix(p_hat, 50),
        main = "Bayesian KDE (diagonal bandwidth)",
        xlab = expression(x[1]), ylab = expression(x[2]))
points(x, pch = 20, cex = 0.4)

Online learning

The isotropic model supports sequential updates: after fitting on an initial batch, you can incorporate new observations cheaply without reprocessing the original data.

n_init <- 100
x_init <- x[1:n_init, ]
x_more <- x[(n_init + 1):n, ]

fit_init   <- ep_kde_isotropic(x_init, prior_shape = 1, prior_rate = 1)
fit_update <- ep_kde_online(x_init, x_more,
                             prior_shape = 1, prior_rate = 1,
                             fit_old = fit_init)
fit_all    <- ep_kde_isotropic(x, prior_shape = 1, prior_rate = 1)

cat("Online posterior mean precision :", fit_update$post_shape / fit_update$post_rate, "\n")
#> Online posterior mean precision : 2.323879
cat("Offline (all-at-once) mean prec :", fit_all$post_shape    / fit_all$post_rate,    "\n")
#> Offline (all-at-once) mean prec : 2.316491