关于scipy.stats.wasserstein_distance无权重p=1实现的技术疑问

Great question! Let's break down how SciPy implements the 1-dimensional Wasserstein distance (for p=1, no weights) step by step, and address all your concerns:

First, you’re correct that for 1-dimensional distributions, the Wasserstein distance ( W_1 ) is equivalent to the integral of the absolute difference between the two cumulative distribution functions (CDFs):

( W_1 = \int_{-\infty}^{+\infty} |U(x) - V(x)| dx )
where ( U(x) ) is the empirical CDF of ( u_values ), and ( V(x) ) is the empirical CDF of ( v_values ). This is a key result from optimal transport theory—for 1D, the optimal transport plan simplifies to matching quantiles, leading directly to this integral form.

Let’s walk through the code and connect it to the integral definition:

1. Sorting and Merging Values (Steps 1-2) #

u_sorter = np.argsort(u_values)
v_sorter = np.argsort(v_values)
all_values = np.concatenate((u_values, v_values))
all_values.sort(kind='mergesort')

The goal here is to find all the critical points where either ( U(x) ) or ( V(x) ) changes. These are exactly the unique values from both distributions. By merging and sorting all values, we split the entire real line into non-overlapping intervals: ( [all_values[0], all_values[1]], [all_values[1], all_values[2]], ..., [all_values[-2], all_values[-1]] ).

In each of these intervals, both ( U(x) ) and ( V(x) ) are constant—there are no sample points from either distribution inside the interval, so the empirical CDF doesn’t increase. This lets us compute the integral as a sum over these intervals, rather than a continuous integral.

2. Calculating Interval Lengths (Step 3) #

deltas = np.diff(all_values)

deltas stores the length of each interval: ( deltas[i] = all_values[i+1] - all_values[i] ). This is the "width" we’ll multiply by the constant ( |U-V| ) in each interval.

3. Computing CDF Values for Each Interval (Steps 4-5) #

u_cdf_indices = u_values[u_sorter].searchsorted(all_values[:-1], 'right')
v_cdf_indices = v_values[v_sorter].searchsorted(all_values[:-1], 'right')
v_cdf = v_cdf_indices / v_values.size
u_cdf = u_cdf_indices / u_values.size

Let’s unpack this:

• u_values[u_sorter] is the sorted version of ( u_values ).
• searchsorted(all_values[:-1], 'right') returns, for each left endpoint of an interval (all_values[:-1]), the number of elements in the sorted ( u_values ) that are less than or equal to that endpoint. The 'right' flag ensures we count all elements ≤ x (since it returns the insertion point after existing entries of x).
• Dividing by the total number of elements (u_values.size) gives us the empirical CDF value ( U(x) ) at that left endpoint—and since ( U(x) ) is constant across the interval, this is the value of ( U(x) ) for all x in the interval.

Same logic applies to v_cdf: it’s the constant value of ( V(x) ) across each interval. This directly answers your question (b): u_cdf and v_cdf are exactly the empirical CDF values of the two distributions over each integration interval.

4. Summing to Compute the Integral (Step 6) #

return np.sum(np.multiply(np.abs(u_cdf - v_cdf), deltas))

This answers your question (a): since ( |U(x)-V(x)| ) is constant over each interval, the integral over that interval is just ( |U-V| \times \text{interval length} ). Summing these products across all intervals gives the exact value of the integral ( \int |U-V| dx )—no approximation needed, because we’ve partitioned the line exactly at all points where the CDFs change.

Short answer: No, it’s completely aligned with the Wasserstein distance definition.

The Wasserstein distance measures the distance between two probability distributions, not between two ordered sequences of samples. The empirical CDF of a sample set depends only on the counts of values in each range, not the order of the samples. Sorting the input values doesn’t change the underlying distribution, so the result is valid.

内容的提问来源于stack exchange，提问作者kam