spatstat 包中的加权中位数

Question

spatstat 包中的 weighted.median() 函数返回“10.5”，当我通过 10、11 和 12 的平均加权分数时。我期待“11”的响应（这是 @ 的输出987654322@ 和 matrixStats::weightedMedian())。

加权中位数的概念对我来说不是很自然。是输出不正确，还是我误解了函数的目的？

x <- c(10, 11, 12)
w <- c( 1,  1,  1)

spatstat::weighted.median(x, w)
#> [1] 10.5
spatstat::weighted.quantile(x, w, probs = .5)
#>  50% 
#> 10.5


matrixStats::weightedMedian(x, w)
#> [1] 11
median(x)
#> [1] 11

^{由reprex package (v0.3.0) 于 2020 年 2 月 23 日创建}

Answer 1

我认为这是包装中的一个缺陷，我会解释原因。

首先，weighted.median 实际上只是调用weighted.quantile 并将probs 向量设置为0.5。但是如果你用你的数据调用weighted.quantile，你会得到非常奇怪的结果：

weighted.quantile(x, w)
#>    0%   25%   50%   75%  100% 
#> 10.00 10.00 10.50 11.25 12.00 

这是不对的。

如果您使用body(weighted.quantile) 查看此函数的主体，并遵循逻辑，则在第 10 行将权重归一化为名为@9​​87654329@ 的变量的方式似乎存在问题。为了正常工作，归一化权重应该是一个与x 长度相同的向量，但从 0 开始到 1 结束，其间的间距与权重成正比。

但如果你看看这是如何计算的：

body(weighted.quantile)[[10]]
#> Fx <- cumsum(w)/sum(w)

你可以看到它不是从 0 开始的。在你的例子中，第一个元素是 0.3333。

为了说明情况，让我们用正确的表达式覆盖这一行。 （首先我们需要解锁绑定才能访问该函数）

unlockBinding("weighted.quantile", asNamespace("spatstat"))
body(weighted.quantile)[[10]] <- substitute(Fx <- (cumsum(w) - min(w))/(sum(w) - min(w)))

现在我们得到了加权分位数的正确结果（包括正确的中位数）

weighted.quantile(x, w)
#>   0%  25%  50%  75% 100% 
#> 10.0 10.5 11.0 11.5 12.0 

Answer 2

这里有一个更基本的问题是关于在小有限样本中分位数的定义（包括中位数）。

R 基函数quantile.default 的帮助文件说有一个参数type，有7 个不同的选项，会给出不同的答案。这些在帮助文件中引用的 Rob Hyndman 的一篇精美文章中进行了详细描述。 quantile.default 的默认值为 type=7。

spatstat::weighted.quantile 中的算法执行type=4 的模拟（根据其帮助文件）；即对累积分布函数F(x)进行线性插值，然后计算逆函数。 该算法在 spatstat 代码中正确实现。

您提到的另一个包中的加权中位数正在计算加权中位数的不同定义。

非常感谢您引起我们对这个示例的关注。这可能会促使我们扩展spatstat::weighted.median 的实现以包含其他类型。

顺便说一句，CRAN 包的错误报告应发布在包的错误报告页面上，如 CRAN 上所示。看到这篇文章真是太幸运了。但是非常感谢你们两位发现这个问题。