使用 numpy 离散接近高斯分布

Question

我正在尝试为 n >= 2 获得与高斯分布的离散接近度。

假设 n = 2，则离散接近度为 [0.5, 0.5]。

当 n = 3 时，它会是 [0.25, 0.5, 0.25]

当 n = 4 时，它会是 [0.125, 0.375, 0.375, 0.125]

希望你明白我的意思。

返回的离散邻近数组总和应始终为 1，因为所有分布。

这是我的代码：

import numpy as np
import matplotlib.pyplot as plt
import math
import scipy 
from random import randint

def discrete_gauss(n):
    g = [0.5, 0.5]
    f = g
    for x in range(1, n - 1):
        f = np.convolve(f,g)

    if(sum(f) != 1):
        print("The distribution sum is not 1.")
    else:
        return f

现在 'discrete_gauss' 在我使用 (1 < n < 68) 时效果很好，但是当我输入 (n > 67) 时，f 的总和不同于 1（有时多有时少），我不知道不知道为什么。有人知道吗？

对于我试图保持简短的凌乱问题，我们深表歉意。如果事情不清楚，我很乐意澄清。谢谢

Answer 1

tl;博士

阅读此 paper 关于使用浮点数学的挑战，然后重新考虑您的方法。

解决方案

这是生成所需 "distribution" 的替代过程，可避免 np.convolve 执行的求和中的浮点舍入错误：

import numpy as np
import scipy.special as sps

def discrete_gauss(n):
    f = np.array([sps.comb(n - 1, i, exact=True) for i in range(n)], dtype='O')
    f = np.float64(f)/np.float64(f).sum()

    if not np.allclose(f.sum(), 1.0):
        raise ValueError("The distribution sum is not close to 1.\n" 
                         "f.sum(): %s" % f.sum())

    return f

解法说明

你要的序列相当于帕斯卡三角的第n层（见Wiki on Binomial theorem上方的图），归一化后可以表示一个概率。上面的解决方案使用标准 Python int 值（在 Python 3 中默认为任意精度）找到第 n 级别的值，然后切换到浮点数学仅在归一化步骤的最后（即 np.float64(f)/np.float64(f).sum()）。

注意在上面的检查中使用 not np.allclose(f.sum(), 1.0)，而不是 f.sum() != 1.0。正如下面 更深入的研究 部分所讨论的，f.sum() 将等于 1.0 对于 n 从 1-1000 的约 90% 的值.但是，一般来说，您不能假设浮点计算的结果与使用实数进行等效计算的结果完全匹配（有关所有详细信息，请参见 paper）。在处理浮点数时，您通常（我的意思是几乎总是）检查结果是否接近（ie 等于给定的 tolerance/error）到您的预期值，而不是不等于它。

深入探讨

这个解决方案并不完美。 n 的大多数值产生的结果总和正好等于 1.0，但有些则不然。以下代码检查 discrete_gauss(n) 的结果以获取 1-1000 的 n 值：

nnot1 = []
for n in range(1,1001):
    if discrete_gauss(n).sum() != 1.0:
        nnot1.append(n)

print('discrete_gauss(n).sum() was not equal to 1.0 for %d values of n.' % len(nnot1))
print(nnot1)

输出：

discrete_gauss(n).sum() was not equal to 1.0 for 75 values of n.
[78, 89, 110, 114, 125, 127, 180, 182, 201, 206, 235, 248, 273, 342, 346, 348, 365, 373, 383, 390, 402, 403, 421, 427, 429, 451, 454, 471, 502, 531, 540, 556, 558, 574, 579, 584, 587, 595, 600, 609, 617, 631, 633, 647, 648, 651, 657, 669, 674, 703, 705, 728, 731, 763, 765, 772, 778, 783, 798, 816, 837, 852, 858, 860, 861, 867, 874, 877, 906, 912, 941, 947, 959, 964, 972]

因此，对于这些值的 ~8%，dicrete_gauss(n).sum() 并不完全等于 1.0。但是，由于没有出现错误，np.allclose(dicrete_gauss(n).sum(), 1.0) 总是 True.

备注

scipy.speical.comb(n, k, exact=True)给出第(n, k)个二项式系数作为标准Pythonint，相当于k中的第k个值14=]Pascal 三角形的第 1 层。

使用 numpy 离散接近高斯分布

discrete proximity to gaussian distribution using numpy

python

numpy

gaussian

convolution

tl;博士

解决方案

解法说明

深入探讨

备注