获取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'来使它工作,至少我这样做...如果我只使用'alpha'和'我会收到错误beta'与Python 2.7.1和SciPy 0.10.0).但即使在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匹配上层分位数的界限大于或等于该分位数中的所有值.