2

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

在此处输入图像描述

其中F^{-1}u是边缘和的累积分布的逆概率分布函数v，源自称为x和的真实数据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，而 Scipyppf可以用来计算逆F^{-1}的？

4

2 回答 2

2

请注意，当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，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)
于 2020-12-07T03:00:26.713 回答
0

我想我有点晚了，但这是我会为一个精确的解决方案做的（只使用 numpy）：

import numpy as np
from numpy.random import randn
n = 100
m = 80
p = 2
x = np.sort(randn(n))
y = np.sort(randn(m))
a = np.ones(n)/n
b = np.ones(m)/m
# cdfs
ca = np.cumsum(a)
cb = np.cumsum(b)

# points on which we need to evaluate the quantile functions
cba = np.sort(np.hstack([ca, cb]))
# weights for integral
h = np.diff(np.hstack([0, cba]))

# construction of first quantile function
bins = ca + 1e-10 # small tolerance to avoid rounding errors and enforce right continuity
index_qx = np.digitize(cba, bins, right=True)    # right=True becouse quantile function is 
                                                 # right continuous
qx = x[index_qx] # quantile funciton F^{-1}      

# construction of second quantile function 
bins = cb + 1e-10 
index_qy = np.digitize(cba, bins, right=True)    # right=True becouse quantile function is 
                                                 # right continuous
qy = y[index_qy] # quantile funciton G^{-1}

ot_cost = np.sum((qx - qy)**p * h)
print(ot_cost)

如果您有兴趣，您可以在这里找到更详细的基于 numpy 的 ot 问题实现，以及双重和原始解决方案：https ://github.com/gnies/1d-optimal-transport 。（虽然我仍在努力）。

于 2021-10-24T17:39:32.553 回答