Python中类似于Matlab的经验性cdf函数

Question

我有一些matlab代码，希望能够重写成python。这是一个简单的程序，计算某个分布并将其以双对数刻度绘制出来。

我的问题出在计算累积分布函数上。以下是matlab代码：

for D = 1:10
    delta = D / 10;
    for k = 1:n
        N_delta = poissrnd(delta^-alpha,1);
        Y_k_delta = ( (1 - randn(N_delta)) / (delta.^alpha) ).^(-1/alpha);
        Y_k_delta = Y_k_delta(Y_k_delta > delta);
        X(k) = sum(Y_k_delta);
        %disp(X(k))

    end
    [f,x] = ecdf(X);

    plot(log(x), log(1-f))
    hold on
end

在Matlab中，我可以简单地使用：

[f,x] = ecdf(X);

要在点x处获取cdf(f)。 这里是相关文档。
在Python中，它更加复杂：

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

alpha = 1.5
n = 1000
X = []
for delta in range(1,5):
    delta = delta/10.0
    for k in range(1,n + 1):
        N_delta = np.random.poisson(delta**(-alpha), 1)
        Y_k_delta = ( (1 - np.random.random(N_delta)) / (delta**alpha) )**(-1/alpha)
        Y_k_delta = [i for i in Y_k_delta if i > delta]
        X.append(np.sum(Y_k_delta))

    ecdf = ECDF(X)

    x = np.linspace(min(X), max(X))
    f = ecdf(x)
    plt.plot(np.log(f), np.log(1-f))

plt.show()

我的图看起来很奇怪，绝对不如matlab的顺畅。
我认为问题在于我不理解 ECDF 函数或它的工作方式与matlab不同。
我使用了这个解决方案（最多点数的），但看起来它并没有正确地工作。

- Photon Light

---

1. 我建议在这些情况下将图形添加到您的问题中，这会很有帮助。
2. 您熟悉Matlab的loglog和plt.loglog吗？

- Andras Deak -- Слава Україні

1. 是的，我知道，我需要一些时间学习如何附加它们...
2. 不，我从未使用过它们。

- Photon Light

1. 只需在编辑器菜单中点击小型景观图标，您就可以通过计算机或链接从网络中包含图像。
2. 我建议您查看这些：它们的工作方式类似于相应的“plot”，但具有更多功能，并且会为您处理“log”。通常情况下，它们甚至更好：它们保留了对数比例尺，同时保留了您的原始变量：如果您的“x”从1到1000，则最终的“xticklabel”将不会写入“log（1000）”，而是“1000”，并且刻度将以对数方式选择。这与绘制“log” vs “log”不同，因此取决于您的需求。

- Andras Deak -- Слава Україні

Answer 1

一旦您获得了样本，您可以使用np.unique*和np.cumsum的组合轻松计算ECDF：

import numpy as np

def ecdf(sample):

    # convert sample to a numpy array, if it isn't already
    sample = np.atleast_1d(sample)

    # find the unique values and their corresponding counts
    quantiles, counts = np.unique(sample, return_counts=True)

    # take the cumulative sum of the counts and divide by the sample size to
    # get the cumulative probabilities between 0 and 1
    cumprob = np.cumsum(counts).astype(np.double) / sample.size

    return quantiles, cumprob

例如：

from scipy import stats
from matplotlib import pyplot as plt

# a normal distribution with a mean of 0 and standard deviation of 1
n = stats.norm(loc=0, scale=1)

# draw some random samples from it
sample = n.rvs(100)

# compute the ECDF of the samples
qe, pe = ecdf(sample)

# evaluate the theoretical CDF over the same range
q = np.linspace(qe[0], qe[-1], 1000)
p = n.cdf(q)

# plot
fig, ax = plt.subplots(1, 1)
ax.hold(True)
ax.plot(q, p, '-k', lw=2, label='Theoretical CDF')
ax.plot(qe, pe, '-r', lw=2, label='Empirical CDF')
ax.set_xlabel('Quantile')
ax.set_ylabel('Cumulative probability')
ax.legend(fancybox=True, loc='right')

plt.show()

---

* 如果您使用的numpy版本低于1.9.0，则np.unique将不接受return_counts关键字参数，您将收到TypeError：

TypeError: unique() got an unexpected keyword argument 'return_counts'

在这种情况下，解决方法是获取“反向”索引集，并使用np.bincount计算出现次数：

quantiles, idx = np.unique(sample, return_inverse=True)
counts = np.bincount(idx)

- ali_m

---

看起来很不错，谢谢！我只是把“return_counts=True”改成了“True”，因为在这个版本中我遇到了错误。顺便说一句，祝你的论文好运！ - Photon Light

如果你遇到了错误 unexpected keyword argument 'return_counts'，那么你可能正在运行一个不支持 np.unique 的 return_counts 参数的较旧版本的 numpy（这是在 v1.9 中添加的）。在这种情况下，调用 np.unique(..., True) 将对应于 np.unique(..., return_index=True)，它将给出唯一值的索引而不是它们的计数作为第二个返回变量，因此你的 CDF 将不正确。请参见我的更新以获取解决方法。 - ali_m