Wasserstein Median of Histograms

Compute a 1D Wasserstein geometric median of midpoint-discrete histograms with common breaks, using a modified Weiszfeld algorithm on exact finite quantile intervals. Return its projection onto the original histogram bins and retain the unprojected finite estimate separately.

Usage

histmed(hists, weights = NULL, L = 2000L, ...)

Arguments

hists

a length-\(N\) list of histograms ("histogram" objects) with identical breaks.

weights

a weight for each histogram; if NULL (default), uniform weights are used. Otherwise, it should be a length-\(N\) vector of nonnegative weights.

L

retained for compatibility and validated as a positive integer. The finite quantile backend uses exact mass intervals, so this argument no longer controls an approximation grid.

...

extra parameters including

abstol

absolute solver tolerance (default: 1e-8).

reltol

relative solver tolerance (default: 1e-8).

maxiter

maximum number of iterations (default: 496).

print.progress

logical; whether to show current iteration (default: FALSE).

Value

a "histogram" display of the median, with the unprojected estimate, solver diagnostics, and projection metadata.

Details

Inputs are discrete measures at bin midpoints, using normalized counts (or the canonical mass field of an estimated histogram). The unprojected finite result is saved in unprojected_measure. For the legacy histogram display it is binned onto the original breaks, using right-closed bins with the lowest endpoint included, and represented at bin midpoints. This projection generally changes the measure; it is not an exact unrestricted Wasserstein summary or geodesic. Downstream histogram functions consistently use the displayed midpoint measure, not the separate unprojected result. Break values must agree exactly across inputs; no grid snapping is done. When an estimated histogram is reused, its previous computational metadata is discarded. Returned diagnostics describe the current calculation. Returned mass and counts are the same probabilities; density equals mass divided by bin width. Counts are not rounded observations, and effective_sample_size is NA because an estimated measure has no implied sample size. The objective/history and solver status refer to the unprojected summary.

Examples

# \donttest{
#----------------------------------------------------------------------
#                      Binned from Two Gaussians
#
# Generate 12 histograms from N(-4,1/4) and 8 from N(4,1/4)
#----------------------------------------------------------------------
# COMMON SETTING
set.seed(100)
bk = seq(from=-10, to=10, length.out=20)
n_signal = 12

# GENERATE HISTOGRAMS WITH COMMON BREAKS
hist_all = list()
for (i in 1:n_signal){
  hist_all[[i]] = hist(stats::rnorm(200, mean=-4, sd=0.5), breaks=bk)
}












for (j in (n_signal+1):20){
  hist_all[[j]] = hist(stats::rnorm(200, mean=+4, sd=0.5), breaks=bk)
}









# COMPUTE THE BARYCENTER AND THE MEDIAN 
h_bary = histbary(hist_all)
h_med  = histmed(hist_all)

# VISUALIZE
xt   <- round(h_med$mids, 1) 
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
barplot(hist_all[[1]]$density, col=rgb(0,0,1,1/4), 
        ylim=c(0, 0.75), main="Two Types", names.arg=xt)
barplot(hist_all[[20]]$density, col=rgb(1,0,0,1/4), 
        ylim=c(0, 0.75), add=TRUE)
barplot(h_med$density,  names.arg=xt, main="Median", ylim=c(0, 0.75))
barplot(h_bary$density, names.arg=xt, main="Barycenter", ylim=c(0, 0.75))

par(opar)
# }