jsd_kde_nd: n-dimensional JSD via multivariate kernel density estimation

In phontrast: Contrast and Separation Metrics for Phonological Categories

n-dimensional JSD via multivariate kernel density estimation

Description

Computes Jensen-Shannon divergence between two categories in an arbitrary n-dimensional acoustic space using multivariate KDE. The default engine uses the ks package; a faster diagonal-Gaussian engine is available for diagonal bandwidths.

Usage

jsd_kde_nd(
  data,
  features,
  group = "category",
  bw = c("Hpi", "Hscv", "Hpi.diag", "scott.diag"),
  eval_on = c("pooled", "group1", "group2", "pooled_sample"),
  eval_n = NULL,
  eval_seed = NULL,
  engine = c("ks", "fast_diag", "fast_diagonal"),
  chunk_size = 1000L,
  method = c("mc", "legacy"),
  density = c("kde", "mvnorm"),
  mc_n = 10000L,
  loo = TRUE,
  bw_scale = 1
)

Arguments

data	A data frame containing observations from exactly two categories.
features	Character vector of column names giving the acoustic dimensions (e.g., MFCC1..MFCC13, F1/F2/duration).
group	String: name of the column giving the category labels (e.g., "vowel", "segment"). Must have exactly two unique values in data.
bw	Bandwidth selection method. One of "Hpi", "Hscv", "Hpi.diag", or "scott.diag". The first three are passed to ks::Hpi(), ks::Hscv(), or ks::Hpi.diag() for multivariate inputs. "scott.diag" uses a diagonal Scott rule-of-thumb bandwidth matrix. For one-dimensional inputs, these map to stats::bw.SJ(), stats::bw.ucv(), stats::bw.nrd0(), and Scott's rule, respectively, with a robust fallback for constant samples.
eval_on	Where to evaluate the KDEs (method = "legacy" only). "pooled" (default) evaluates on all observations from both categories; "group1" or "group2" evaluate on the respective group only. "pooled_sample" evaluates on a sampled subset of pooled observations and requires eval_n. Ignored when method = "mc" (which always evaluates each category at its own observations).
eval_n	Optional positive integer giving the maximum number of evaluation points to use. If supplied, evaluation points are sampled from the set chosen by eval_on.
eval_seed	Optional integer seed used only when eval_n causes evaluation-point subsampling. If NULL, the current R random-number state is used.
engine	KDE evaluation engine. "ks" uses ks::kde(). "fast_diag" uses a chunked diagonal-Gaussian evaluator and requires bw = "scott.diag" or bw = "Hpi.diag" for multivariate KDE. "fast_diagonal" is accepted as an alias for "fast_diag".
chunk_size	Positive integer controlling the number of evaluation points processed per chunk by engine = "fast_diag".
method	Estimator: "mc" (default) for the Monte-Carlo plug-in estimate of the continuous JSD, or "legacy" for the pre-1.2.0 self-normalized sample-point index. Ignored when density = "mvnorm".
density	Density model behind the estimate: "kde" (default) estimates each category's density by kernel density estimation; "mvnorm" fits one multivariate normal per category and estimates the continuous JSD between the two Gaussians by Monte-Carlo (no closed form exists). Under "mvnorm" the KDE-specific arguments (bw, engine, eval_on, chunk_size, method, eval_n, loo) do not apply; the Monte-Carlo sample size is set by mc_n and eval_seed makes the draw reproducible.
mc_n	Positive integer; number of Monte-Carlo samples drawn from each fitted Gaussian when density = "mvnorm" (default 10000). The estimator draws mc_n fresh points from each category's fitted Gaussian and averages the log density ratio, so it targets the JSD between the two fitted Gaussians rather than a resubstitution estimate at the observed points. Larger values reduce Monte-Carlo variance. Ignored when density = "kde".
loo	Logical; if TRUE (default) the Monte-Carlo estimator uses a partial leave-one-out correction on each category's self-density to reduce resubstitution bias. The correction removes a sample-size-scaled fraction n / (n + 20) of each point's own kernel: half at 20 tokens per category (the min_tokens default), approaching the full leave-one-out correction as the category grows. Removing only part of the self-kernel keeps the corrected density strictly positive at isolated points, so small but real divergences remain small positive values rather than being floored to exactly 0 (as the full leave-one-out correction did through phontrast 2.0.2). Ignored when method = "legacy".
bw_scale	Positive number multiplying the selected kernel bandwidth on the standard-deviation scale: univariate bandwidths are multiplied by bw_scale and bandwidth matrices by bw_scale^2. The default 1 uses the selected bandwidth unchanged; 0.5 and 2 give the halved and doubled bandwidths of the smoothing-sensitivity check in rank_contrasts(). Ignored when density = "mvnorm".

Details

By default (method = "mc") JSD is estimated with a Monte-Carlo plug-in: each category's KDE is evaluated at that category's own observations and the log density ratio against the mixture is averaged. This is a consistent estimator of the continuous JSD in any dimension. method = "legacy" reproduces the pre-1.2.0 self-normalized sample-point estimate (a bounded relative separation index rather than the continuous JSD); use it only to reproduce results from phonJSD 1.0.0.

Value

A single numeric JSD value in bits, bounded in [0, 1].

Examples

set.seed(2026)
vowels <- data.frame(
  vowel = rep(c("ih", "eh"), each = 40),
  f1 = c(rnorm(40, 500, 55), rnorm(40, 565, 60)),
  f2 = c(rnorm(40, 1980, 150), rnorm(40, 1870, 155))
)

# One-dimensional JSD, for example a single formant or duration.
jsd_kde_nd(vowels, features = "f1", group = "vowel")

# Two-dimensional JSD in F1/F2 space.
jsd_kde_nd(vowels, features = c("f1", "f2"), group = "vowel")

# Faster high-dimensional path: diagonal Scott bandwidth and sampled
# pooled evaluation points.
jsd_kde_nd(
  vowels,
  features = c("f1", "f2"),
  group = "vowel",
  bw = "scott.diag",
  eval_n = 40,
  eval_seed = 2026,
  engine = "fast_diag"
)

phontrast documentation built on Oct. 7, 2026, 5:06 p.m.