Python: оценка плотности ядра для положительных значений

Question

Я хочу получить оценку плотности ядра для положительных точек данных. Используя пакет Python Scipy Stats, я разработал следующий код.

def get_pdf(data):
    a = np.array(data)
    ag = st.gaussian_kde(a)
    x = np.linspace(0, max(data), max(data))
    y = ag(x)
    return x, y

Это прекрасно работает для большинства наборов данных, но дает ошибочный результат для "всех положительных" точек данных. Чтобы убедиться, что это работает правильно, я использую численное интегрирование для вычисления площади под этой кривой.

def trapezoidal_2(ag, a, b, n):
    h = np.float(b - a) / n
    s = 0.0
    s += ag(a)[0]/2.0
    for i in range(1, n):
        s += ag(a + i*h)[0]
    s += ag(b)[0]/2.0
    return s * h

Так как данные распространяются в области (0, int (max (data))), мы должны получить значение, близкое к 1, при выполнении следующей строки.

b = 1
data = st.pareto.rvs(b, size=10000)
data = list(data)

a = np.array(data)
ag = st.gaussian_kde(a)
trapezoidal_2(ag, 0, int(max(data)), int(max(data))*2)

Но при тестировании он дает значение, близкое к 0,5.

Но когда я интегрирую от -100 до максимума (данные), он дает значение, близкое к 1.

trapezoidal_2(ag, -100, int(max(data)), int(max(data))*2+200)

Причина в том, что ag (KDE) определено для значений меньше 0, даже если исходный набор данных содержит только положительные значения.

Итак, как я могу получить оценку плотности ядра, которая учитывает только положительные значения, так что область под кривой в области (o, max (данные)) близка к 1?

Answer 1

Выбор пропускной способности очень важен при оценке плотности ядра. Я думаю, что правило Скотта и правило Сильвермана хорошо работают для распределения, подобного гауссовскому. Тем не менее, они не работают хорошо для распределения Парето.

Цитата из документа :

Выбор полосы пропускания сильно влияет на оценку, полученную из KDE (намного больше, чем фактическая форма ядра). Выбор пропускной способности может быть сделано по эмпирическому правилу, путем перекрестной проверки, с помощью плагина методы "или другими способами; см. [3] , [4] для обзоров. gaussian_kde использует эмпирическое правило, по умолчанию используется правило Скотта.

Попробуйте использовать другие значения пропускной способности, например:

import numpy as np
import matplotlib.pyplot as plt

from scipy import stats

b = 1

sample = stats.pareto.rvs(b, size=3000)
kde_sample_scott = stats.gaussian_kde(sample, bw_method='scott')
kde_sample_scalar = stats.gaussian_kde(sample, bw_method=1e-3)


# Compute the integrale:
print('integrale scott:', kde_sample_scott.integrate_box_1d(0, np.inf))
print('integrale scalar:', kde_sample_scalar.integrate_box_1d(0, np.inf))

# Graph:
x_span = np.logspace(-2, 1, 550)
plt.plot(x_span, stats.pareto.pdf(x_span, b), label='theoretical pdf')
plt.plot(x_span, kde_sample_scott(x_span), label="estimated pdf 'scott'")
plt.plot(x_span, kde_sample_scalar(x_span), label="estimated pdf 'scalar'")
plt.xlabel('X'); plt.legend();

дает:

integrale scott: 0.5572130540733236
integrale scalar: 0.9999999999968957

и

Мы видим, что kde, использующий метод Скотта, неверен.