Python 中的一维 Wasserstein 距离

下面的公式是 Wasserstein 距离/最优传输的特例，当源和目标分布为x and y（也称为边际分布）是一维的，即向量。

where F^{-1}是边际累积分布的逆概率分布函数u and v，源自称为x and y，均由正态分布生成：

import numpy as np
from numpy.random import randn
import scipy.stats as ss

n = 100
x = randn(n)
y = randn(n)

公式中的积分如何用python和scipy编码？我猜测 x 和 y 必须转换为排名边际，它们是非负的并且总和为 1，而 Scipy 的ppf可以用来计算逆F^{-1}'s?

请注意，当n变大我们有一组排序的n样本接近以 1/n、2/n、...、n/n 采样的逆 CDF。例如。：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
plt.plot(norm.ppf(np.linspace(0, 1, 1000)), label="invcdf")
plt.plot(np.sort(np.random.normal(size=1000)), label="sortsample")
plt.legend()
plt.show()

另请注意，从 0 到 1 的积分可以近似为 1/n、2/n、...、n/n 的总和。

因此我们可以简单地回答你的问题：

def W(p, u, v):
    assert len(u) == len(v)
    return np.mean(np.abs(np.sort(u) - np.sort(v))**p)**(1/p)

请注意，如果len(u) != len(v)您仍然可以应用线性插值方法：

def W(p, u, v):
    u = np.sort(u)
    v = np.sort(v)
    if len(u) != len(v):
        if len(u) > len(v): u, v = v, u
        us = np.linspace(0, 1, len(u))
        vs = np.linspace(0, 1, len(v))
        u = np.linalg.interp(u, us, vs)
    return np.mean(np.abs(u - v)**p)**(1/p)

如果您有关于数据分布类型的先验信息，但没有其参数，则另一种方法是找到数据的最佳拟合分布（例如，使用scipy.stats.norm.fit） 对彼此而言u and v然后以所需的精度进行积分。例如。：

from scipy.stats import norm as gauss
def W_gauss(p, u, v, num_steps):
    ud = gauss(*gauss.fit(u))
    vd = gauss(*gauss.fit(v))
    z = np.linspace(0, 1, num_steps, endpoint=False) + 1/(2*num_steps)
    return np.mean(np.abs(ud.ppf(z) - vd.ppf(z))**p)**(1/p)