python - 在 Python 和 NumPy 中量化正态分布 float

Question

让数组A中的值从高斯采样分配。我想将 A 中的每个值替换为 n_R 之一R 中的“representatives”，以便总量化误差为最小化。

这是进行线性量化的 NumPy 代码:

n_A, n_R = 1_000_000, 256
mu, sig = 500, 250
A = np.random.normal(mu, sig, size = n_A)
lo, hi = np.min(A), np.max(A)
R = np.linspace(lo, hi, n_R)
I = np.round((A - lo) * (n_R - 1) / (hi - lo)).astype(np.uint32)

L = np.mean(np.abs(A - R[I]))
print('Linear loss:', L)
-> Linspace loss: 2.3303939600700603

虽然这有效，但量化误差很大。有没有更聪明的怎么办？我认为人们可以利用 A正态分布或可能使用迭代过程最小化“损失”函数。

更新在研究这个问题时，我发现了 related question关于“加权”量化。调整他们的方法有时会给出更好的量化结果:

from scipy.stats import norm

dist = norm(loc = mu, scale = sig)
bounds = dist.cdf([mu - 3*sig, mu + 3*sig])
pp = np.linspace(*bounds, n_R)
R = dist.ppf(pp)

# Find closest matches
lhits = np.clip(np.searchsorted(R, A, 'left'), 0, n_R - 1)
rhits = np.clip(np.searchsorted(R, A, 'right') - 1, 0, n_R - 1)

ldiff = R[lhits] - A
rdiff = A - R[rhits]
I = lhits
idx = np.where(rdiff < ldiff)[0]
I[idx] = rhits[idx]

L = np.mean(np.abs(A - R[I]))
print('Gaussian loss:', L)
-> Gaussian loss: 1.6521974945326285

K-means 聚类可能更好，但似乎太慢了适用于大型阵列。

Answer 1

K-均值

K-means clustering might be better but seem to be too slow to be practical on large arrays.

对于一维聚类情况，有比 K 均值更快的算法。请参阅https://stats.stackexchange.com/questions/40454/determine-different-clusters-of-1d-data-from-database

我选择了其中一种算法，Jenks Natural Breaks ，并在数据集的随机子样本上运行它:

A_samp = np.random.choice(A, size=10000)
breaks = np.array(jenkspy.jenks_breaks(A_samp, n_classes=n_R))
R = (breaks[:-1] + breaks[1:]) / 2

这相当快，整个数据集的量化损失约为 1.28。

为了可视化这些方法的作用，我绘制了每个方法针对断点 R 内的索引得出的断点的 cdf。

根据定义，高斯是一条直线。这意味着它在分布的每个百分位数上都有相同数量的中断。线性方法在分布的中间花费很少的中断，而在尾部使用大部分中断。 Jenks 在两人之间找到了妥协方案。

自动搜索更低的损失

看着上面的图表，我有了一个想法:所有这些选择中断的方法在分位数域中绘制时都是各种 S 型曲线。 (如果你把它看作一个真正延伸的 sigmoid，那么高斯就很适合。)

我编写了一个函数，该函数使用单个变量“强度”对每条曲线进行参数化，强度是 sigmoid 弯曲的速度。完成后，我使用 scipy.optimize.minimize 自动搜索一条使损失最小化的曲线。

事实证明，如果让 Scipy 对此进行优化，它会选择一条非常接近 Jenks 的曲线强度，并且它找到的曲线比 Jenks 稍差，损失约为 1.33。

您可以看到采用此失败方法的笔记本 here .

使用 2^16 浮点进行量化

在需要创建 2^16 个不同代表的情况下，使用 Jenks 在计算上是不可行的。但是，您可以做一些非常接近的事情:具有少量类的 Jenks 加上线性插值。

这是此操作的代码:

import itertools


def pairwise(iterable):
    "s -> (s0, s1), (s1, s2), (s2, s3), ..."
    a, b = itertools.tee(iterable)
    next(b, None)
    return zip(a, b)


def linspace_jenks(A, n_R, jenks_classes, dist_lo, dist_hi):
    assert n_R % jenks_classes == 0, "jenks_classes must be divisor of n_R"
    simplify_factor = n_R // jenks_classes
    assert jenks_classes ** 2 <= len(A), "Need more data to estimate"
    breaks = jenkspy.jenks_breaks(A, n_classes=jenks_classes)
    # Adjust lowest and highest break to match highest/lowest observed value
    breaks[0] = dist_lo
    breaks[-1] = dist_hi
    linspace_classes = []
    for lo, hi in pairwise(breaks):
        linspace_classes.append(np.linspace(lo, hi, simplify_factor, endpoint=False))
    linspace_classes = np.hstack(linspace_classes)
    assert len(linspace_classes) == n_R
    return linspace_classes

调用示例:

A_samp = np.random.choice(A, size = 2**16)
jenks_R = linspace_jenks(A_samp, n_R, 128, np.min(A), np.max(A))

与线性方法相比，性能如何？在我的系统上，n_R=2^16 的线性损失为 0.009421。下图显示了 linspace_jenks 方法针对 jenks_classes 的每个值获得的损失。

只需 32 个 Jenks 类，并用线性插值填充其余部分，损失就会降至 0.005031。