numpy/scipy 等效于 R ecdf(x)(x) 函数？答案

【问题标题】：numpy/scipy equivalent of R ecdf(x)(x) function?numpy/scipy 等效于 R ecdf(x)(x) 函数？
【发布时间】：2013-03-25 10:13:43
【问题描述】：

在 numpy 或 scipy 中，R 的 ecdf(x)(x) 函数在 Python 中的等价物是什么？是ecdf(x)(x)基本一样吗：

import numpy as np
def ecdf(x):
  # normalize X to sum to 1
  x = x / np.sum(x)
  return np.cumsum(x)

还是需要其他东西？

编辑如何控制ecdf 使用的垃圾箱数量？

【问题讨论】：

This 应该会有所帮助。

标签： python r numpy scipy

【解决方案1】：

试试这些链接：

statsmodels.ECDF

ECDF in python without step function?

示例代码

import numpy as np
from statsmodels.distributions.empirical_distribution import ECDF
import matplotlib.pyplot as plt

data = np.random.normal(0,5, size=2000)

ecdf = ECDF(data)
plt.plot(ecdf.x,ecdf.y)

【讨论】：

scipy 无法做到这一点吗？

【解决方案2】：

这位作者有一个非常好的用户编写的 ECDF 函数示例：John Stachurski's Python lectures。他的系列讲座面向计算经济学研究生；但是，对于任何学习 Python 通用科学计算的人来说，它们都是我的首选资源。

编辑：这已经一岁了，但我想我仍然会回答您问题的“编辑”部分，以防您（或其他人）仍然觉得它有用。

ECDF 确实没有像直方图那样的“bins”。如果 G 是使用数据向量 Z 形成的经验分布函数，则 G(x) 实际上是 Z

有趣的奖励：如果您需要从非常大的流数据创建一个类似于 ECDF 的小尺寸对象，您应该查看 McDermott 等人的这篇“Data Skeletons”论文。

【讨论】：

链接现在已断开。作者将他的python讲座移至：quant-econ.net/py/index.html
新链接也坏了。也许这是新的（希望是永久的）链接？ python.quantecon.org/intro.html

【解决方案3】：

ecdf 的 OP 实现是错误的，你不应该 cumsum() 值。所以不是ys = np.cumsum(x)/np.sum(x)，而是ys = np.cumsum(1 for _ in x)/float(len(x))或者更好的ys = np.arange(1, len(x)+1)/float(len(x))

如果你对额外的依赖没有问题，你可以选择statmodels 的ECDF，或者提供你自己的实现。见下文：

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.distributions.empirical_distribution import ECDF
%matplotlib inline

grades = (93.5,93,60.8,94.5,82,87.5,91.5,99.5,86,93.5,92.5,78,76,69,94.5,
          89.5,92.8,78,65.5,98,98.5,92.3,95.5,76,91,95,61)


def ecdf_wrong(x):
    xs = np.sort(x) # need to be sorted
    ys = np.cumsum(xs)/np.sum(xs) # normalize so sum == 1
    return (xs,ys)
def ecdf(x):
    xs = np.sort(x)
    ys = np.arange(1, len(xs)+1)/float(len(xs))
    return xs, ys

xs, ys = ecdf_wrong(grades)
plt.plot(xs, ys, label="wrong cumsum")
xs, ys = ecdf(grades)
plt.plot(xs, ys, label="handwritten", marker=">", markerfacecolor='none')
cdf = ECDF(grades)
plt.plot(cdf.x, cdf.y, label="statmodels", marker="<", markerfacecolor='none')
plt.legend()
plt.show()

【讨论】：

【解决方案4】：

R 中的ecdf 函数返回经验累积分布函数，因此完全等价的应该是：

def ecdf(x):
    x = np.sort(x)
    n = len(x)
    def _ecdf(v):
        # side='right' because we want Pr(x <= v)
        return (np.searchsorted(x, v, side='right') + 1) / n
    return _ecdf

np.random.seed(42)
X = np.random.normal(size=10_000)
Fn = ecdf(X)
Fn([3, 2, 1]) - Fn([-3, -2, -1])
## array([0.9972, 0.9533, 0.682 ])

如图所示，它给出了正态分布的correct 68–95–99.7% probabilities。

【讨论】：