假设数组A
中的值是从高斯分布中采样的。我想用R
中的n_R
“代表”之一替换A
中的每个值,以便使总量化误差最小化。
下面是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聚类可能更好,但似乎太慢,无法在大型数组上实用。
2条答案
按热度按时间2guxujil1#
K-means聚类可能更好,但似乎太慢,无法在大型数组上实用。
对于1D聚类的情况,有比K均值更快的算法。参见https://stats.stackexchange.com/questions/40454/determine-different-clusters-of-1d-data-from-database
我选择了其中一个算法Jenks Natural Breaks,并在数据集的随机子样本上运行它:
字符串
这是相当快的,并且对于完整数据集得到约1.28的量化损失。
为了直观地显示这些方法的作用,我绘制了每个方法产生的break的cdf与break的R内的索引。
x1c 0d1x的数据
根据定义,高斯是一条直线。这意味着它在分布的每个百分位数处具有相等数量的中断。线性方法在分布的中间花费很少的中断,并且在尾部使用大部分中断。詹克斯在他们两人之间找到了一个折中的办法。
自动寻找损失更小的
看着上面的图表,我有了一个想法:当在分位数域中绘制时,所有这些选择中断的方法都是各种S形曲线。(如果你把它看作是一个真正伸展的S形,那么高斯模型就很适合。)
我写了一个函数,用一个变量强度来参数化每一条曲线,强度是S形曲线应该弯曲的速度。一旦我有了它,我就使用
scipy.optimize.minimize
自动搜索一条最小化损失的曲线。事实证明,如果你让Scipy优化它,它会选择一个非常接近詹克斯的曲线强度,它发现的曲线比Jenks的曲线略差,损失约1.33。
你可以在这里看到这个失败方法的笔记本。
sbtkgmzw2#
部分是为了新奇,主要是为了完整性,我演示了@Homer512正确建议的可能性-MILP实现。我期望它的准确性是优秀的,它的表现是介于“差”和“可怕”之间。
我用一个非常小的问题来演示,这样当您调试和查看约束矩阵时,它们是清晰的,并且我的RAM不会爆炸。
字符串