使用 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 npimport matplotlib.pyplot as pltimport mathimport scipy from random import randintdef 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现在，当我使用 (1 < n < 68) 时，'discrete_gauss' 效果很好，但是当我输入 (n > 67) 时，f 的总和与 1 不同（有时多有时少），我不知道为什么。任何人都有任何线索？对不起，我试图保持简短的凌乱问题。如果事情不清楚，我很乐意澄清。谢谢。

Answer 1

阅读这篇关于使用浮点数学的挑战的论文，然后重新考虑你的方法。

解决方案

这是生成所需“分布”的替代过程，可避免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帕斯卡三角形的第 th 层（参见二项式定理 Wiki顶部的图），归一化以便它可以表示概率。上述解决方案使用标准 Pythonint值（在 Python 3 中默认为任意精度）来查找n第 th 级中的值，然后仅在归一化步骤（即np.float64(f)/np.float64(f).sum()）的最后切换到浮点数学。

请注意在not np.allclose(f.sum(), 1.0)上面的检查中使用，而不是f.sum() != 1.0。正如下面更深入的潜水部分所讨论的f.sum()，1.0对于n1-1000的值的约 90%将等于。但是，通常您不能假设浮点计算的结果与使用实数进行等效计算得到的结果完全匹配（有关所有详细信息，请参阅本文）。在处理浮点数时，您通常（我的意思是几乎总是）检查结果是否接近（即等于在给定的容差/误差范围内）您的预期值，而不是等于它。

更深的潜水

这个解决方案并不完美。大多数值n产生的结果正好等于1.0，但有些则不是。以下代码检查1-1000 之间的discrete_gauss(n)值的结果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)th 二项式系数作为标准 Python int，它相当于帕斯卡三角形k的n第 th 层中的th 值。