дискретная близость к гауссову распределению с использованием 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; dr

Прочтите эту статью о трудностях использования математики с плавающей запятой, а затем пересмотрите свой подход.

Решение

Вотальтернативная процедура для генерации желаемого «распределения», позволяющая избежать ошибок округления с плавающей запятой при суммировании, которое выполняет 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 -го уровня треугольника Паскаля (см. рисунок вверху Вики по теореме о биноме ), нормализованного таким образом, чтобы он мог представлять вероятность.Вышеупомянутое решение использует стандартные значения 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 для ~ 90% значений n от 1-1000.Однако, как правило, вы не можете предполагать, что результат вычисления с плавающей запятой будет точно совпадать с результатом, который вы получите при эквивалентном вычислении с использованием действительных чисел (см. Эту бумагу для всех подробностей кровавых).При работе с числами с плавающей точкой обычно (под этим я имею в виду почти всегда) проверяйте, что результат близок (, т.е. равен заданному допуску / ошибке) к ожидаемому значению, а не равен ему.

Более глубокое погружение

Это решение не идеально.Большинство значений 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)-й биномиальный коэффициенткак стандартный Python int, который эквивалентен k th значению в n th уровне треугольника Паскаля.