NumPy 百分位函数不同于 MATLAB 的百分位函数

Posted 2023-03-12

【中文标题】NumPy 百分位函数不同于 MATLAB 的百分位函数

【英文标题】：NumPy percentile function different from MATLAB's percentile function

【发布时间】：2014-09-06 01:08:54

【问题描述】：

当我尝试在 MATLAB 中计算第 75 个百分位数时，我得到的值与在 NumPy 中不同。

MATLAB：

>> x = [ 11.308 ;   7.2896;   7.548 ;  11.325 ;   5.7822;   9.6343;
     7.7117;   7.3341;  10.398 ;   6.9675;  10.607 ;  13.125 ;
     7.819 ;   8.649 ;   8.3106;  12.129 ;  12.406 ;  10.935 ;
    12.544 ;   8.177 ]

>> prctile(x, 75)

ans =

11.3165

Python + NumPy：

>>> import numpy as np

>>> x = np.array([ 11.308 ,   7.2896,   7.548 ,  11.325 ,   5.7822,   9.6343,
     7.7117,   7.3341,  10.398 ,   6.9675,  10.607 ,  13.125 ,
     7.819 ,   8.649 ,   8.3106,  12.129 ,  12.406 ,  10.935 ,
    12.544 ,   8.177 ])

>>> np.percentile(x, 75)
11.312249999999999

我也用 R 检查了答案，我得到了 NumPy 的答案。

R：

> x <- c(11.308 ,   7.2896,   7.548 ,  11.325 ,   5.7822,   9.6343,
+          7.7117,   7.3341,  10.398 ,   6.9675,  10.607 ,  13.125 ,
+          7.819 ,   8.649 ,   8.3106,  12.129 ,  12.406 ,  10.935 ,
+         12.544 ,   8.177)
> quantile(x, 0.75)
     75% 
11.31225 

这里发生了什么？有没有办法让 Python 和 R 的行为反映 MATLAB 的？

【问题讨论】：

您能告诉我们 MATLAB 使用的公式吗？ R 有9 different ways 来计算分位数。似乎 MATLAB 回答了 R 中的第 75 个匹配 quantile(x, 0.75, type=2) 和 quantile(x, 0.75, type=5)。

当然——来自 MATLAB 帮助页面：（我不能在这里评论它，因为它太长了）mathworks.com/help/stats/prctile.html（你可能需要展开底部的“算法”按钮）

【参考方案1】：

MATLAB 显然默认使用中点插值。 NumPy 和 R 默认使用线性插值：

In [182]: np.percentile(x, 75, interpolation='linear')
Out[182]: 11.312249999999999

In [183]: np.percentile(x, 75, interpolation='midpoint')
Out[183]: 11.3165

---

了解linear和midpoint之间的区别，考虑这个简单的例子：

In [187]: np.percentile([0, 100], 75, interpolation='linear')
Out[187]: 75.0

In [188]: np.percentile([0, 100], 75, interpolation='midpoint')
Out[188]: 50.0

---

编译最新版本的 NumPy（使用 Ubuntu）：

mkdir $HOME/src
git clone https://github.com/numpy/numpy.git
git remote add upstream https://github.com/numpy/numpy.git
# Read ~/src/numpy/INSTALL.txt
sudo apt-get install libatlas-base-dev libatlas3gf-base
python setup.py build --fcompiler=gnu95
python setup.py install

使用git 而不是pip 的优点是升级（或降级）到其他版本的 NumPy 非常容易（并且您也可以获得源代码）：

git fetch upstream
git checkout master # or checkout any other version of NumPy
cd ~/src/numpy
/bin/rm -rf build
cdsitepackages    # assuming you are using virtualenv; otherwise cd to your local python sitepackages directory
/bin/rm -rf numpy numpy-*-py2.7.egg-info
cd ~/src/numpy
python setup.py build --fcompiler=gnu95
python setup.py install

【讨论】：

嗯——你使用的是什么版本的 NumPy？我没有添加插值关键字参数的选项。 （我试过 numpy 和 scipy.stats.scoreatpercentile）

@James: The interpolation parameter 是在 1.9.0 版本中添加的。

pip install -U numpy 应该可以工作 AFAIK。但是，我发布了一种使用 git 编译/安装最新版本的替代方法。

边缘情况并非如此。 Matlab 对分配给最低和最高数据点的百分位数之间的所有内容使用中点插值。对于超出此范围的百分位数，它会分配数据集的最小值和最大值。 source

这不是不完整的，它是错误。 Matlab 默认使用线性插值（见nl.mathworks.com/help/stats/prctile.html）而不是中点。它只是与 numpy 不同的线性插值变体。有关线性插值的不同变体，请参阅 en.wikipedia.org/wiki/Percentile 和 en.wikipedia.org/wiki/Quantile。

【参考方案2】：

由于即使在@cpaulik 发表评论后接受的答案仍然不完整，所以我在这里发布希望是更完整的答案（尽管出于简洁的原因，并不完美，见下文）。

使用 np.percentile(x, p, interpolation='midpoint') 只会对非常具体的值给出相同的答案，即当 p/100 是 1/n 的倍数时，n 是元素的数量的数组。在最初的问题中，确实如此，因为 n=20 和 p=75，但总的来说这两个函数不同。

Matlab 的 prctile 函数的简短仿真如下：

def quantile(x,q):
    n = len(x)
    y = np.sort(x)
    return(np.interp(q, np.linspace(1/(2*n), (2*n-1)/(2*n), n), y))

def prctile(x,p):
    return(quantile(x,np.array(p)/100))

这个函数，作为 Matlab 的函数，给出了一个从 min(x) 到 max(x) 的分段线性输出。 Numpy 的百分位函数，插值=中点，返回两个最小元素的平均值和两个最大元素的平均值之间的分段常量函数。在原始问题中为数组绘制两个函数会给出the picture in this link （抱歉无法嵌入）。红色虚线标记了 75% 的百分位，这两个函数实际上重合。

附：这个函数实际上不等同于 Matlab 的原因是它只接受一维 x，对更高维的东西给出错误。另一方面，Matlab 接受更高的 dim x 并在第一个（非平凡的）维度上运行，但正确实现它可能需要更长的时间。然而，这个函数和 Matlab 的函数都应该正确地处理 p / q 的更高维输入（感谢使用 np.interp 来处理它）。

【讨论】：

嗨，Marco，我对您的回答很感兴趣，但是通过您的功能，我得到了静态输出。

对不起，静态输出是什么意思？此函数接受一维数组和一个数字作为输入并输出一个数字，例如：prctile(np.arange(100),1) 输出 0.5。它不适用于高维数组。如果您将多维 q （或类似数组）传递给它，它将返回一个相同形状的数组，例如prctile(np.arange(100),[1, 3]) 输出 [0.5, 2.5]。