python - R 和 Python 中 `medcouple()` 实现的区别

Question

我发现R 和Python 中的medcouple() 实现之间存在差异。考虑一个由 10 组成的数组，重复 480 次并以 [1,2,3,4,5,6,7,8,9,10,11,12,13,14 结尾,15,16,17,18,19,20,21]。对于这个数组，R 和 Python 中的 medcouple() 返回不同的结果。

以下 R 代码返回 0:

library(mrfDepth)
values = c(rep(10, 480), 
c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21))
medcouple(values, FALSE)

但是下面的Python代码:

from statsmodels.stats.stattools import medcouple
arr=[10.0]*480 + [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]
medcouple(arr)
returns 1!

恕我直言 R 就在这里，但是有没有原创开发者可以对此发表评论？

Answer 1

查看您的数据，可以预期 medcouple 很小，但呈阳性。 Medcouple 是偏度的量度。正如Hubert and Vandervieren所述:

It clearly follows from this definition that the medcouple always lies between −1 and 1. A distribution that is skewed to the right has a positive value for the medcouple, whereas the MC becomes negative at a left skewed distribution. Finally, a symmetric distribution has a zero medcouple.

在您的数据中，您有很多 10，加上左侧的 9 个值(1 到 9)和右侧的 11 个值(11 到 21)。所以，它一点向右倾斜。

在你的计算中，由于四舍五入，它们返回零(我检查了你的 python 代码，它返回 0，而不是 1。)但是如果你在数据中放的 10 不是那么多，你可以看到一个小的正值:

> medcouple([10.0]*3 + [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21])
array(0.05263158)

更新

实际上，我最初的解释是错误的......对此感到抱歉!

我检查了您提供的 R 代码，它返回 1，而不是您看到的零。我还绘制了(使用 R)不同样本的 medcouple(向初始数据添加不同数量的 10s)。

我不仅使用了 mrfDepth 库，还使用了同样实现 medcouple 统计的 robustbase。他们返回了相同的值。

为了帮助理解这种不对称度量的行为，我添加了一些使用偏度统计的图。

最后，为了让事情更有趣，我将所有这些结果与一个非常相似的数据样本进行了比较(只省略了 2 个数字:20 和 21)。

见下面的代码:

library(robustbase) # mc (also a medcouple implementation)
library(mrfDepth) # medcouple
library(moments) # skewness

symmetric.sample = function(n) {
    c(c(1:19), rep(10,n))
}
skewed.sample = function(n) {
    c(c(1:21), rep(10,n))
}

xlab = "# 10s added"
n = 1:150

png("skew.png", width=600, height=1000, pointsize=8, res=160)
par(mfrow=c(4,2), pch=20)
hist(sapply(30, symmetric.sample), breaks=0:21, xlab="symmetric sample (with 30 10s added)", main="")
hist(sapply(30, skewed.sample), breaks=0:21, xlab="skewed sample (with 30 10s added)", main="")
plot(n, sapply(sapply(n, symmetric.sample), robustbase::mc), col="red", xlab=xlab, ylab="robustbase's mc on symmetric sample")
plot(n, sapply(sapply(n, skewed.sample), robustbase::mc), col="red", xlab=xlab, ylab="robustbase's mc on skewed sample")
plot(n, sapply(sapply(n, symmetric.sample), medcouple, do.reflect=FALSE), col="red", xlab=xlab, ylab="mrfDepth's mc on symmetric sample")
plot(n, sapply(sapply(n, skewed.sample), medcouple, do.reflect=FALSE), col="red", xlab=xlab, ylab="mrfDepth's mc on skewed sample")
plot(n, sapply(sapply(n, symmetric.sample), skewness), col="red", xlab=xlab, ylab="skewness on symmetric sample")
plot(n, sapply(sapply(n, skewed.sample), skewness), col="red", xlab=xlab, ylab="skewness on skewed sample")
dev.off()

现在，关于 python 的 stattools 结果，结果大不相同。添加 4 个或更多 10 做原始数据(1 到 21)，medcouple 返回 0。

我测试了Jordi Gutiérrez Hermoso's python implementation .它与 R medcouples 函数一致。