使用 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 67) 时，f 的总和不同于 1（有时更多有时更少），我不知道为什么. 有人知道吗？

抱歉，我试图简短地回答这个凌乱的问题。如果不清楚，我会很乐意澄清。 谢谢。

Answer 1

tl;博士

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

解决方案

这是生成所需“分布”的替代过程，可避免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

解法说明

你想要的序列等价于帕斯卡三角形的nth 级（参见Wiki on Binomial theorem 顶部的图），经过归一化处理后可以表示概率。上述解决方案使用标准 Python int 值（在 Python 3 中默认为任意精度）来查找 nth 级别中的值，然后仅在归一化步骤的最后切换到浮点数学（即np.float64(f)/np.float64(f).sum())。

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

深潜

这个解决方案并不完美。 n 的大多数值产生的结果与1.0 的总和完全一致，但有些则不然。下面的代码检查 discrete_gauss(n) 的结果中 n 从 1 到 1000 的值：

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)th 二项式系数作为标准 Python int 给出，这相当于帕斯卡三角形的 nth 级别中的 kth 值。