как интегрировать распределение Гаусса с массивом дисперсии

Question

Я пытаюсь вычислить интеграцию распределения Гаусса.У меня есть массив сигма.

sigma = np.array([0.2549833 , 0., 0.42156247, 0. , 0.,  0., 0.79124217, 0.54235005, 0.79124217, 0.        ,     0.        , 0.32532629, 0.46753655, 0.60605513, 0.55420338,      0.        , 0.38053264, 0.42690288, 0.        , 0.63526099])

И формула распределения Гаусса:

def gaussian(x, mu, sig):
if sig != 0:
    return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))

Интеграция этого распределения Гаусса:

I = np.zeros(len(sigma), dtype=float)
for i in range(0, len(sigma)):
    I[i] = quad(gaussian(x, mu = 0, sig = sigma[i]), 0, 105)

, но это не работаетпотому что quad функция выдает ошибку.Как я могу получить массив Интеграции в этом случае?

Answer 1

Интеграция гауссовского pdf дает кумулятивную функцию плотности, cdf.Итак, что вам нужно, это scipy.stats.norm.cdf.И это поддерживает векторизованные оценки, поэтому нет необходимости в цикле for в вашем коде.

Обратите внимание, что cdf интегрируется из отрицательной бесконечности.Кроме того, параметр сигма называется scale.

Answer 2

Проблема была в использовании вами quad.Вот правильная версия.Несколько замечаний:

• Вы создаете I длины, равной количеству значений сигмы, но вы ничего не делали, когда сигма равна 0. Итак, теперь я воссоздал sigmaиметь только ненулевые значения.Таким образом, я удалил условие if sig != 0:

• Правильное использование в соответствии с документами заключается в том, что вы сначала передаете определение функции, затем пределы, а затем - аргументы функции.используя ключевое слово args.

• quad вернет кортеж, содержащий интеграл и абсолютную ошибку.Следовательно, вы должны использовать индекс [0], чтобы получить целое значение для хранения в I.Вы можете сохранить абсолютную ошибку в другом списке, если хотите.

Для вашего удобства я добавил график значений Integral против sigma.

Модифицированная часть вашего кода

sigma = sigma[sigma !=0]

def gaussian(x, mu, sig):
        return np.exp(-(x - mu)**2/ (2 * sig**2))

I = np.zeros(len(sigma), dtype=float)

for i in range(0, len(sigma)):
    I[i]  = quad(gaussian, 0, 105,  args=(0, sigma[i]))[0]