Теоретическая функция нормального распределения в scipy

Question

Мне нужно построить нормальное кумулятивное распределение для заданных ребер бинов:

bin_edges = np.array([1.02,  4.98,  8.93, 12.89, 16.84, 20.79, 24.75, 28.7])
mean = 15.425
standard_deviation = 6.159900567379315

Сначала я сделал:

cdf = ((1 / (np.sqrt(2 * np.pi) * standard_deviation)) *
   np.exp(-0.5 * (1 / standard_deviation * (bin_edges - mean))**2))
cdf = cdf.cumsum()
cdf /= cdf[-1]

Другой способ, который я нашел:

cdf = scipy.stats.norm.cdf(bin_edges, loc=mean, scale=standard_deviation)

Вывод этих двух методов должен быть одинаковым, но это не так:

First: [0.0168047  0.07815162 0.22646339 0.46391741 0.71568769 0.89247475 
0.97468339 1.]
Second: [0.0096921  0.04493372 0.14591031 0.34010566 0.59087116 0.80832701
0.93495018 0.98444529]

Для меня это выглядит как результат scipy cdf () хуже.Что я делаю не так?

Answer 1

Проблема

Вы пытаетесь вычислить CDF на каждом ребре ячейки, вычисляя значение следующего интеграла на каждом ребре ячейки:

Причина, по которой ваш результат не совпадает с результатом scipy, заключается в том, что scipy делает интеграцию лучше, чем вы.Вы эффективно интегрируете нормальный PDF путем суммирования по области «столбцов» гистограммы, которые эффективно определяют ваши bin_edges.Это не даст достаточно точный результат, пока количество бинов не станет намного, намного выше (возможно, по крайней мере, в тысячах).Ваш подход к нормализации также не подходит, поскольку на самом деле вам нужно делить на интеграл PDF от -inf до inf, а не от 1.02 до 28.7.

С другой стороны,Numpy просто вычисляет высокоточное численное приближение решения интеграла в замкнутой форме.Используемая функция называется scipy.special.ndtr.Вот это реализация в коде Scipy .

Решение

Вместо интегрирования путем суммирования площадей столбцов вы можете выполнить фактическое численное интегрирование от -inf до x, чтобы получить результат с точностью, приближающейся к точности scipy.stats.norm.cdf,Вот код для того, как это сделать:

import scipy.integrate as snt

def pdf(x, mean, std):
    return ((1/((2*np.pi)**.5 * std)) * np.exp(-.5*((x - mean)/std)**2))

cdf = [snt.quad(pdf, -np.inf, x, args=(mean, std))[0] for x in bin_edges]

Версия Scipy ndtr написана на C, но вот примерное приближение Python для сравнения:

import scipy.special as sps

def ndtr(x, mean, std):
    return .5 + .5*sps.erf((x - mean)/(std * 2**.5))

Тестирование

import scipy.special as sps
import scipy.stats as sts
import scipy.integrate as snt

bin_edges = np.array([1.02,  4.98,  8.93, 12.89, 16.84, 20.79, 24.75, 28.7])
mean = 15.425
std = 6.159900567379315

with np.printoptions(linewidth=9999):
    print(np.array([snt.quad(pdf, -np.inf, x, args=(mean, std))[0] for x in bin_edges]))
    print(ndtr(bin_edges, mean, std))
    print(sts.norm.cdf(bin_edges, loc=mean, scale=std))

Вывод:

[0.00968036 0.04497664 0.14584988 0.34034101 0.59084202 0.80811081 0.93496465 0.98442171]
[0.00968036 0.04497664 0.14584988 0.34034101 0.59084202 0.80811081 0.93496465 0.98442171]
[0.00968036 0.04497664 0.14584988 0.34034101 0.59084202 0.80811081 0.93496465 0.98442171]

Таким образом, при точной интеграции результаты с помощью метода, который вы использовали, с высокой точностью соответствуют результатам scipy.stats.norm.cdf.