Python: Как получить интегральную функцию распределения для непрерывных значений данных?

Question

У меня есть набор значений данных, и я хочу получить CDF (накопительную функцию распределения) для этого набора данных.

Поскольку это непрерывная переменная, мы не можем использовать подход биннинга, как упомянуто в ( Как правильно получить интегральную функцию распределения для моих данных в python? ). Поэтому я придумал следующий подход.

import scipy.stats as st

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

def get_cdf(data):
    a = np.array(data)
    ag = st.gaussian_kde(a)

    cdf = [0]
    x = []
    k = 0

    max_data = max(data)

    while (k < max_data):
        x.append(k)
        k = k + 1

    sum_integral = 0

    for i in range(1, len(x)):
        sum_integral = sum_integral + (trapezoidal_2(ag, x[i - 1], x[i], 2))
        cdf.append(sum_integral)

    return x, cdf

Вот как я использую этот метод.

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

В идеале я должен получить значение, близкое к 1, в конце списка y_cdf. Но я получаю значение, близкое к 0,57.

Что здесь не так? Мой подход правильный?

Спасибо.

Answer 1

Я думаю, что это просто:

def get_cdf(data):
  return sorted(data), np.linspace(0, 1, len(data))

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

Когда я сравниваю это с аналитическим результатом, я получаю то же самое:

x_cdf, y_cdf = get_cdf(st.pareto.rvs(1, size=10000))

import matplotlib.pyplot as plt
plt.semilogx(x_cdf, y_cdf)
plt.semilogx(x_cdf, st.pareto.cdf(x_cdf, 1))

Answer 2

Значение cdf в x является интегралом pdf между -inf и x, но вы вычисляете его между 0 и x. Возможно, вы предполагаете, что pdf равен 0 для x <0, но это не так: </p>

rs = np.random.RandomState(seed=52221829)
b = 1
data = st.pareto.rvs(b, size=10000, random_state=rs)
ag = st.gaussian_kde(data)

x = np.linspace(-100, 100)
plt.plot(x, ag.pdf(x))

Так что, вероятно, здесь происходит что-то не так: вы не проверяете свои предположения.

Ваш код для вычисления интеграла мучительно медленный, есть лучшие способы сделать это с помощью scipy, но gaussian_kde предоставляет метод integrate_box_1d для интеграции pdf. Если вы берете интеграл из -inf, все выглядит правильно.

cdf = np.vectorize(lambda x: ag.integrate_box_1d(-np.inf, x))
plt.plot(x, cdf(x))

Интегрируя от 0 до x, вы получаете то же, что вы видите сейчас (справа от 0), но это совсем не cdf:

wrong_cdf = np.vectorize(lambda x: ag.integrate_box_1d(0, x))
plt.plot(x, wrong_cdf(x))

Answer 3

Не уверен, почему ваша функция не работает точно, но один из способов вычисления CDF заключается в следующем:

def get_cdf_1(data):

    # start with sorted list of data
    x = [i for i in sorted(data)]

    cdf = []

    for xs in x:
        # get the sum of the values less than each data point and store that value
        # this is normalised by the sum of all values
        cum_val = sum([i for i in data if i <= xs])/sum(data) 
        cdf.append(cum_val)

    return x, cdf

Нет никаких сомнений в том, что более быстрый способ вычислить это, используя массивы, а не добавление значений ксписок, но он возвращает значения в том же формате, что и исходный пример.