Distribution representations and transport objectives

One observation is a probability distribution

Use ot_measure() for finite weighted atoms. A numeric vector gives scalar locations; a matrix gives one multivariate location per row. Probability masses within an observation and weights across observations answer different questions.

mu <- ot_measure(c(0, 2), mass = c(1, 3))
nu <- as_ot_measure(ecdf(c(0, 0, 4)))
data.frame(location = mu$support[, 1], probability = mu$mass)
#>   location probability
#> 1        0        0.25
#> 2        2        0.75
observations <- ot_collection(list(first = mu, second = nu), weights = c(2, 1))
observations
#> Collection of 2 probability measures in 1 dimension(s)
#>   label outer_weight atoms
#>   first    0.6666667     2
#>  second    0.3333333     2

Here mu places probabilities 1/4 and 3/4 at its two locations. The collection gives the whole distribution mu weight 2/3. Duplicating every observation in an empirical sample leaves that empirical probability distribution unchanged; it does not automatically increase its outer weight. Input masses and outer weights must be finite, nonnegative, and have positive total. Zero masses and duplicate support rows retain their original indices in a requested plan.

A cost and a distance have different units

For order p, the Euclidean transport cost is distance raised to p, and the Wasserstein distance is the pth root of the optimal cost. With p = 2, support in meters gives a cost in squared meters and a distance in meters.

comparison <- ot_distance(mu, nu, p = 2, return_plan = TRUE)
comparison
#> Transport comparison | 1d | euclidean | order 2 
#> Status: SUCCESS ( finite_monotone_transport )
#> Transport cost: 3 
#> Wasserstein distance: 1.732051
stopifnot(comparison$status == "success")
c(row_error = max(abs(rowSums(comparison$plan) - mu$mass)),
  column_error = max(abs(colSums(comparison$plan) - nu$mass)))
#>    row_error column_error 
#> 0.000000e+00 5.551115e-17

Automatic method selection uses monotone finite transport for scalar supports and network simplex otherwise. A numerical answer is an exact optimization target evaluated in floating-point arithmetic. Inspect status and residuals before treating a computed plan as optimal.

A supplied cross-matrix must declare whether it contains distances or already powered costs. Its rows and columns refer to the original source and target indices. A metric declaration is supplied by the caller; the package cannot verify a triangle inequality from one cross-matrix alone.

D <- abs(outer(mu$support[, 1], nu$support[, 1], "-"))
from_distance <- ot_distance(mu, nu, p = 2,
                            cost = ot_cost(D, type = "distance", metric = TRUE))
from_cost <- ot_distance(mu, nu, p = 2,
                        cost = ot_cost(D^2, type = "cost", metric = TRUE))
stopifnot(isTRUE(all.equal(from_distance$transport_cost, from_cost$transport_cost)))

For arbitrary nonmetric costs, keep metric = FALSE. The result then reports a transport cost without claiming that it defines a Wasserstein metric.

Entropy regularization changes the objective

With plan entries P, the Sinkhorn target is

\[ \langle C,P\rangle + \varepsilon\sum_{ij} P_{ij}(\log P_{ij}-1), \]

where a zero entry contributes zero. epsilon has the same units as the powered cost. Scaling coordinates by a positive factor s requires scaling epsilon by s^p to preserve the corresponding regularized problem.

regularized <- ot_distance(mu, nu, method = "sinkhorn", epsilon = 0.5)
regularized
#> Transport comparison | sinkhorn | euclidean | order 2 
#> Status: SUCCESS ( marginal_tolerance )
#> Transport cost: 3 
#> Rooted transport component: 1.732051 
#> Regularized objective: 1.961222
stopifnot(regularized$status == "success")
regularized[c("transport_cost", "transport_cost_root", "regularized_objective")]
#> $transport_cost
#> [1] 3
#> 
#> $transport_cost_root
#> [1] 1.732051
#> 
#> $regularized_objective
#> [1] 1.961222

The transport component is evaluated at the regularized plan. Its root is not claimed to be an exact Wasserstein distance or a metric. The full regularized objective can be negative, so taking its square root is not a median construction. return_plan = FALSE reduces the returned object, but internal dense allocations may still be needed.

Scalar barycenters and medians

The barycenter minimizes sum(weights * W2^2). The median minimizes sum(weights * W2). Both use inner Wasserstein order two. The second problem is a geometric median of complete quantile functions, not separate pointwise quantile medians.

barycenter <- ot_summary(observations, target = "barycenter")
median <- ot_summary(observations, target = "median")
barycenter
#> Wasserstein barycenter | unrestricted_1d | quantile 
#> Status: SUCCESS ( quantile_average )
#> Objective: 0.6666667 | estimated atoms: 3
median
#> Wasserstein median | unrestricted_1d | quantile 
#> Status: SUCCESS ( collision_optimality )
#> Objective: 0.5773503 | estimated atoms: 2
stopifnot(barycenter$status == "success", median$status == "success")
data.frame(location = barycenter$estimate$support[, 1],
           mass = barycenter$estimate$mass)
#>   location      mass
#> 1 0.000000 0.2500000
#> 2 1.333333 0.4166667
#> 3 2.666667 0.3333333

The finite engine combines all cumulative-mass breakpoints and integrates the resulting quantile intervals exactly, up to floating-point arithmetic. It needs neither a uniform probability grid nor simulated repeated samples. The median uses modified Weiszfeld updates and verifies a collision before accepting an input as a solution. An unrepresentable positive mass interval produces an explicit error rather than silently disappearing.

Choose a conversion deliberately

Input	Conversion	Meaning
Samples or weighted atoms	ot_measure()	Finite atomic probability measure
Base ECDF or fitted weighted CDF	as_ot_measure()	Exact finite jumps, preserving ties
Histogram	as_ot_measure()	Normalized bin counts at bin midpoints
Density values on a grid	ot_density()	Trapezoidal node masses, a finite approximation
Image	ot_image()	Normalized pixel masses at explicit pixel centers
Gaussian parameters	ot_gaussian()	A separate parametric Gaussian object

A histogram density height is not a probability mass for unequal-width bins. The input converter uses counts; displayed fitted densities equal mass divided by bin width. Histogram summary displays can differ from the separately retained unprojected finite estimate.

grid <- seq(-3, 3, length.out = 31)
density_measure <- ot_density(grid, dnorm(grid))
image_measure <- ot_image(matrix(c(1, 0, 2, 1, 0, 0), 2, 3, byrow = TRUE))
density_measure$metadata$approximation
#> [1] "finite discretization of supplied density curve"
image_measure$metadata[c("x", "y", "order")]
#> $x
#> [1] 0.5 1.5 2.5
#> 
#> $y
#> [1] 1.5 0.5
#> 
#> $order
#> [1] "row-major"

The density calculation normalizes the discretized curve on the supplied grid; truncation and grid refinement remain modeling choices. Image intensities are pixel masses, even with nonuniform coordinates. Normalizing discards total intensity. Default image coordinates increase rightwards and upwards, with row-major indexing. Legacy image functions use other documented coordinate conventions; see the migration tutorial.

Work directly with Gaussian parameters

g1 <- ot_gaussian(mean = 0, covariance = 1)
g2 <- ot_gaussian(mean = 4, covariance = 9)
ot_gaussian_distance(g1, g2)
#> Transport comparison | gaussian_formula | euclidean | order 2 
#> Status: SUCCESS ( gaussian_parameter_formula )
#> Transport cost: 20 
#> Wasserstein distance: 4.472136
gaussbary1d(c(0, 4), c(1, 9), weights = c(1, 3))[c("mean", "var", "objective", "status")]
#> $mean
#> [1] 3
#> 
#> $var
#> [1] 6.25
#> 
#> $objective
#> [1] 3.75
#> 
#> $status
#> [1] "success"
gaussmed1d(c(0, 1, 2), c(4, 16, 9))[c("mean", "var", "objective", "status")]
#> $mean
#> [1] 1.211325
#> 
#> $var
#> [1] 10.31261
#> 
#> $objective
#> [1] 1.115355
#> 
#> $status
#> [1] "success"

Univariate Gaussian W2 geometry is Euclidean geometry of (mean, standard deviation). Consequently the barycenter averages these two parameters and squares the averaged standard deviation to recover variance. The median is their joint geometric median. Zero univariate variance represents a point mass.

For multivariate summaries, gaussbarypd() and gaussmedpd() take an n by p mean matrix and a p by p by n covariance array. Every covariance must be symmetric and strictly positive definite, including when using these pd functions with p = 1. No ridge or PSD projection is applied. General covariance medians use reweighted Gaussian barycenter subproblems and report numerical stopping criteria without a global-optimality claim. These parameter objects are not automatically converted into random finite samples.