Barycenter of Histograms by Benamou et al. (2015)

Given multiple histograms represented as "histogram" S3 objects, compute Wasserstein barycenter. We need one requirement that all histograms in an input list hists must have same breaks. See the example on how to construct a histogram on predefined breaks/bins. Break values must agree exactly across inputs; no grid snapping is done.

Usage

histbary15B(hists, p = 2, weights = NULL, lambda = NULL, ...)

Arguments

hists

a length-\(N\) list of histograms ("histogram" object) of same breaks.

p

an exponent for the order of the distance (default: 2).

weights

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

lambda

a regularization parameter; if NULL (default), a paper's suggestion would be taken, or it should be a strictly positive real number.

...

extra parameters including

abstol: stopping criterion for iterations (default: 1e-8).
init.vec: a strictly positive initial mass vector (default: uniform mass).
maxiter: maximum number of iterations (default: 496).
nthread: compatibility control; only 1 is supported.
print.progress: a logical to show current iteration (default: FALSE).

Value

a "histogram" object of barycenter.

Details

Each input is the discrete measure at its bin midpoints with normalized bin counts, including unequal-width bins. The entropy coefficient lambda multiplies the sum of transport-plan entries times their logarithms; the cost is midpoint distance raised to p. The finite regularization coefficient changes the objective from the unregularized barycenter target. The shared solver now uses log-domain Sinkhorn gradients and objective-checked mirror descent; historical names remain compatibility entry points. Returned mass and counts are probability masses, and density is mass divided by bin width. No sample-size-dependent rounding is performed; effective_sample_size is NA. Status, accepted objective history, and inner diagnostics are available in fit. The full objective uses log(plan) - 1; this differs from log(plan) by a constant and has the same minimizer.

References

Benamou J, Carlier G, Cuturi M, Nenna L, Peyré G (2015). “Iterative Bregman Projections for Regularized Transportation Problems.” SIAM Journal on Scientific Computing, 37(2), A1111-A1138. ISSN 1064-8275, 1095-7197, doi:10.1137/141000439 .

Examples

#----------------------------------------------------------------------
#                      Binned from Two Gaussians
#
# EXAMPLE : Very Small Example for CRAN; just showing how to use it!
#----------------------------------------------------------------------
# GENERATE FROM TWO GAUSSIANS WITH DIFFERENT MEANS
set.seed(100)
x  = stats::rnorm(1000, mean=-4, sd=0.5)
y  = stats::rnorm(1000, mean=+4, sd=0.5)
bk = seq(from=-10, to=10, length.out=20)

# HISTOGRAMS WITH COMMON BREAKS
histxy = list()
histxy[[1]] = hist(x, breaks=bk, plot=FALSE)
histxy[[2]] = hist(y, breaks=bk, plot=FALSE)

# COMPUTE
hh = histbary15B(histxy, maxiter=5)
#> Warning: histbary15B finished with status 'iteration_limit': maximum_iterations. Inspect result$fit.

# VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2))
barplot(histxy[[1]]$density, col=rgb(0,0,1,1/4), 
        ylim=c(0, 0.75), main="Two Histograms")
barplot(histxy[[2]]$density, col=rgb(1,0,0,1/4), 
        ylim=c(0, 0.75), add=TRUE)
barplot(hh$density, main="Barycenter",
        ylim=c(0, 0.75))

par(opar)