Python：如何获得连续数据值的累积分布函数？

Question

我有一组数据值，我想获取该数据集的 CDF（累积分布函数）。

由于这是一个连续变量，我们不能使用 (How to get cumulative distribution function correctly for my data in 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)

理想情况下，我应该在 y_cdf 列表的末尾得到一个接近 1 的值。但我得到的值接近 0.57。

这里出了什么问题？我的方法正确吗？

谢谢。

Answer 1

x 处的 cdf 值是 -inf 和 x 之间 pdf 的积分，但您是在 0 和 x 之间计算它。也许您假设 x

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 2

不确定为什么您的函数无法正常工作，但计算 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

毫无疑问，使用 numpy 数组而不是将值附加到列表中是一种更快的计算方法，但这会以与原始示例相同的格式返回值。

Answer 3

我认为只是：

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))