获取 SciPy 分位数以匹配 Stata xtile 函数

Question

我继承了一些旧的 Stata 代码 (Stata11)，它使用 xtile 函数按分位数对向量中的观察值进行分类（在这种情况下，只是标准的 5 个五分位数，20%、40%、60%、 80%、100%）。

我正在尝试在 Python 中复制一段代码，并且我正在使用 SciPy.stats.mstats 函数 mquantiles() 进行计算。

据我从 Stata 文档和在线搜索中得知，Stata xtile 方法尝试反转数据的经验 CDF，并使用 CDF 平坦的所有观察值的等加权平均值来进行切点。这似乎是对分位数进行分类的一种非常糟糕的方法，但事实就是如此，我相信在某些情况下这是正确的做法。

我的问题是如何让mquantiles() 产生同样的打破惯例。我注意到这个函数有两个参数，alphap 和betap（文档称它们为alpha 和beta，但你需要额外的'p'才能让它工作，至少我这样做......我如果我只在 Python 2.7.1 和 SciPy 0.10.0 中使用“alpha”和“beta”，则会出现错误）。但即使在 SciPy 文档中，我也看不出这些参数的组合是否会在平坦的 CDF 范围内产生平均值。

我看到了作为这个范围的中位数或众数计算的选项，但不是平均值（也不清楚这些具有 alpha 和 beta 的 SciPy 中位数/众数选项是否被计算为 观察或产生平坦 CDF 值的范围。）

任何帮助消除这些不同选项的歧义并找到一些帮助我在 Python 中重新创建 Stata 约定的文档都会很棒。请不要回答只是说“编写自己的分位数函数”。首先，这无助于我理解 Stata 或 SciPy 的约定，其次，鉴于这些数值库，编写自己的分位数函数应该是最后的手段。我当然可以，但如果我需要，那就不好了。

Answer 1

scipy.stats.mquantiles 文档很差，并且在某些地方错误，现在已修复，因此可能会有所帮助... http://docs.scipy.org/scipy/docs/scipy.stats.mstats_basic.mquantiles/。当您指出 alpha/beta、alphap/betap 差异时，该过程就开始了。谢谢。

mquantiles 的实现遵循 R。

最大的区别在于 R 有 9 种离散类型，因为 scipy.stats.mquantiles 从 'alphap' 和 'betap' 计算 'm'，所以 scipy 有一个连续范围的“类型”（因为没有更好的词）。

我承认我并不了解所涉及的统计数据的所有细节，因此我决定采用蛮力评估。我在http://www.biostat.sdu.dk/~biostat/StataReferenceManual/StataRef.pdf 找到了一个 xtile 示例，并且能够将结果与 alphap=0.5 和 betap=0.5（分段线性）相匹配。不是确定的，也不是详尽无遗的，但我现在拥有的一切。

In [1]: import scipy.stats as st

In [9]: st.mstats.mquantiles([23,56,67,123,99,17],prob=[0.5],alphap=0.5,betap=.5)
Out[9]: array([ 61.5])

In [10]: st.mstats.mquantiles([23,56,67,123,99,17],prob=[0.33,0.66],alphap=0.5,betap=.5)
Out[10]: array([ 38.84,  81.72])

In [11]: st.mstats.mquantiles([23,56,67,123,99,17],prob=[0.25,0.5,0.75],alphap=0.5,betap=.5)
Out[11]: array([ 23. ,  61.5,  99. ])

最后一个有点问题，因为其中两个分割点正好在数据集中的值上。 Stata/xtile（至少在我发现的示例中）没有给出分位数的分割点，而是给出了分位数本身。给定排序后的数据集 [17,23,56,67,99,123]，Stata/xtile 将分类为 [1,1,2,3,3,4]，这意味着对于 scipy.stat.mquantiles 匹配上分位数的界限大于或等于该分位数中的所有值。