De-heaping density estimation with adheaping
In adheaping: Characteristic-Function De-Heaping Density Estimation

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 3.5)
library(adheaping)

Heaping and the identifiability limit

Heaped and rounded data concentrate on a coarse grid of round numbers and defeat kernel density estimation, which renders the rounding marks as spurious modes. Rounding to a grid of width D is exactly convolution of the density with a width-D box followed by lattice sampling, so in the characteristic-function domain the density is recovered by dividing out the known box (a continuous generalization of Sheppard's correction) inside the grid-Nyquist band |w| < pi/D. Beyond that band the density is not identifiable from the heaped data alone.

A worked example

We draw a bimodal sample, round it to a grid of 0.5, and compare the naive kernel estimate with the tuning-free combined de-heaping estimator.

set.seed(20260627)
n <- 4000
x <- ifelse(runif(n) < 0.5, rnorm(n, -1.2, 0.5), rnorm(n, 1.2, 0.5))
D <- 0.5
y <- D * round(x / D)
grid <- seq(-6, 6, length.out = 2048)

f_true  <- 0.5 * dnorm(grid, -1.2, 0.5) + 0.5 * dnorm(grid, 1.2, 0.5)
f_naive <- naive_kde(y, grid)
f_adk   <- adkde(y, D, grid)
attr(f_adk, "pick")   # which component the band-capacity gate selected

plot(grid, f_true, type = "l", lwd = 2, xlab = "x", ylab = "density", xlim = c(-3, 3))
lines(grid, f_naive, col = "grey55", lty = 3, lwd = 1.6)
lines(grid, as.numeric(f_adk), col = "#b2182b", lwd = 1.8)
legend("topright", c("truth", "naive KDE", "adkde (combined)"),
       col = c("black", "grey55", "#b2182b"), lty = c(1, 3, 1), lwd = 1.8, bty = "n")

The naive estimate carries the rounding comb; the combined estimator recovers the smooth bimodal density.

Reading the grid and the heaped fraction

The grid and the fraction of a sample that is heaped are read directly from the rounding comb, and the comb is treated as a group-matched atom in the spectral basis by a fourth-order detector, which abstains rather than guessing and so fires only in a minority of cases.

heap_grid(y, grid, near = D)                 # locate the grid from the comb tooth
heap_detect(y = D * round(x / D), span = c(-12.8, 12.8))$D_hat  # NA when it abstains

Baselines

Faithful base-R replicas of the measurement-error deconvolution and Heitjan-Rubin multiple-imputation methods are provided (deconv_kde, heitjan_mi), and the real Kernelheaping stochastic EM is wrapped by sem_kde when that package is installed.