R语言 是否存在加权.median()函数?

xdnvmnnf  于 2022-12-20  发布在  其他
关注(0)|答案(9)|浏览(267)

我正在寻找一些类似于weighted.mean()的形式。我已经通过搜索找到了一些解决方案,写出了整个函数,但会欣赏一些更友好的用户。

hmtdttj4

hmtdttj41#

以下软件包都有计算加权中位数的函数:“芳香.光”、“等渗”、“利马”、“cwhmisc”、“ergm”、“laeken”、“矩阵统计”、“PSCBS”和“bigvis”(在github上)。
为了找到它们,我使用了'sos'包中非常宝贵的findFn(),它是R内置帮助的扩展。

findFn('weighted median')

或者,
第一个月
as???是快捷方式,与?some.functionhelp(some.function)的快捷方式相同

zzwlnbp8

zzwlnbp82#

使用@wkmor1和@Jaitropmange的答案的一些经验。
我已经检查了3个软件包中的3个函数,isotonelaekenmatrixStats。只有matrixStats工作正常。其他两个(就像median(rep(x, times=w)解决方案一样)给出整数输出。只要我计算人口的中位年龄,小数位就很重要。

可重现的例子:人口中位年龄的计算

df <- data.frame(age = 0:100,
                 pop = spline(c(4,7,9,8,7,6,4,3,2,1),n = 101)$y)

library(isotone)
library(laeken)
library(matrixStats)

isotone::weighted.median(df$age,df$pop)
# [1] 36
laeken::weightedMedian(df$age,df$pop)
# [1] 36
matrixStats::weightedMedian(df$age,df$pop)
# [1] 36.164
median(rep(df$age, times=df$pop))
# [1] 35

总结

matrixStats::weightedMedian()是可靠的解决方案

nxowjjhe

nxowjjhe3#

使用相同长度的(整数)权重向量w计算向量x的加权中值:

median(rep(x, times=w))
brtdzjyr

brtdzjyr4#

这只是一个简单的解决方案,几乎可以随时随地使用。

weighted.median <- function(x, w) {
  w <- w[order(x)]
  x <- x[order(x)]

  prob <- cumsum(w)/sum(w)
  ps <- which(abs(prob - .5) == min(abs(prob - .5)))
  return(x[ps])
}
bcs8qyzn

bcs8qyzn5#

很老的帖子,但我只是偶然发现它,并做了一些测试的不同方法. spatstat::weighted.median()似乎是大约14倍的速度比median(rep(x, times=w)),它实际上是显而易见的,如果你想运行函数超过几次.测试是与一个相对较大的调查,约15,000人.

uqdfh47h

uqdfh47h6#

也可以使用stats::density创建加权PDF,然后将其转换为CDF,如here所述:

my_wtd_q = function(x, w, prob, n = 4096) 
  with(density(x, weights = w/sum(w), n = n), 
       x[which.max(cumsum(y*(x[2L] - x[1L])) >= prob)])

那么my_wtd_q(x, w, .5)将是加权中位数。
还可以更小心地通过重新归一化来确保density下的总面积为1。

6bc51xsx

6bc51xsx7#

base 中获得 * 加权中位数 * 的一种方法是按值排序,构建权重的cumsum,并获得权重的sum * 0.5处的值。

medianWeighted <- function(x, w, q=.5) {
  n <- length(x)
  i <- order(x)
  w <- cumsum(w[i])
  p <- w[n] * q
  j <- findInterval(p, w)
  Vectorize(function(p,j) if(w[n] <= 0) NA else
    if(j < 1) x[i[1]] else
      if(j == n) x[i[n]] else
        if(w[j] == p) (x[i[j]] + x[i[j+1]]) / 2 else
          x[i[j+1]])(p,j)
}

使用简单的输入数据将得到以下结果。

medianWeighted(c(10, 40), c(1, 2))
#[1] 40
median(rep(c(10, 40), c(1, 2)))
#[1] 40

medianWeighted(c(10, 40), c(2, 1))
#[1] 10
median(rep(c(10, 40), c(2, 1)))
#[1] 10

medianWeighted(c(10, 40), c(1.5, 2))
#[1] 40
medianWeighted(c(10, 40), c(3, 4))
#[1] 40
median(rep(c(10, 40), c(3, 4)))
#[1] 40

medianWeighted(c(10, 40), c(1.5, 1.5))
#[1] 25
medianWeighted(c(10, 40), c(3, 3))
#[1] 25
median(rep(c(10, 40), c(3, 3)))
#[1] 25

medianWeighted(c(10, 40), c(0, 1))
#[1] 40
medianWeighted(c(10, 40), c(1, 0))
#[1] 10
medianWeighted(c(10, 40), c(0, 0))
#[1] NA

它也可用于其它铁路

medianWeighted(1:10, 10:1, seq(0, 1, 0.25))
[1]  1  2  4  6 10

与其他方法进行比较。

#Functions from other Answers
weighted.median <- function(x, w) {
  w <- w[order(x)]
  x <- x[order(x)]

  prob <- cumsum(w)/sum(w)
  ps <- which(abs(prob - .5) == min(abs(prob - .5)))
  return(x[ps])
}

my_wtd_q = function(x, w, prob, n = 4096) 
  with(density(x, weights = w/sum(w), n = n), 
       x[which.max(cumsum(y*(x[2L] - x[1L])) >= prob)])

weighted.quantile <- function(x, w, probs = seq(0, 1, 0.25),
                              na.rm = FALSE, names = TRUE) {

  if (any(probs > 1) | any(probs < 0)) stop("'probs' outside [0,1]")

  if (length(w) == 1) w <- rep(w, length(x))
  if (length(w) != length(x)) stop("w must have length 1 or be as long as x")

  if (isTRUE(na.rm)) {
    w <- x[!is.na(x)]
    x <- x[!is.na(x)]
  }

  w <- w[order(x)] / sum(w)
  x <- x[order(x)]

  cum_w <- cumsum(w) - w * (1 - (seq_along(w) - 1) / (length(w) - 1))
  res <- approx(x = cum_w, y = x, xout = probs)$y

  if (isTRUE(names)) {
    res <- setNames(res, paste0(format(100 * probs, digits = 7), "%"))
  }
  res
}

方法

M <- alist(
  medRep = median(rep(DF$x, DF$w)),
 isotone = isotone::weighted.median(DF$x, DF$w),
 laeken = laeken::weightedMedian(DF$x, DF$w),
 spatstat1 = spatstat.geom::weighted.median(DF$x, DF$w, type=1),
 spatstat2 = spatstat.geom::weighted.median(DF$x, DF$w, type=2),
 spatstat4 = spatstat.geom::weighted.median(DF$x, DF$w, type=4),
 survey = survey::svyquantile(~x, survey::svydesign(id=~1, weights=~w, data=DF), 0.5)$x[1],
 RAndres = weighted.median(DF$x, DF$w),
 matrixStats = matrixStats::weightedMedian(DF$x, DF$w),
 MichaelChirico = my_wtd_q(DF$x, DF$w, .5),
 Leonardo = weighted.quantile(DF$x, DF$w, .5),
 GKi = medianWeighted(DF$x, DF$w)
)

结果

DF <- data.frame(x=c(10, 40), w=c(1, 2))
sapply(M, eval)
#        medRep        isotone         laeken      spatstat1      spatstat2 
#      40.00000       40.00000       40.00000       40.00000       25.00000 
#     spatstat4         survey        RAndres    matrixStats MichaelChirico 
#      17.50000       40.00000       10.00000       30.00000       34.15005 
#  Leonardo.50%            GKi 
#      25.00000       40.00000 

DF <- data.frame(x=c(10, 40), w=c(1, 1))
sapply(M, eval)
#        medRep        isotone         laeken      spatstat1      spatstat2 
#      25.00000       25.00000       40.00000       10.00000       10.00000 
#     spatstat4         survey        RAndres    matrixStats MichaelChirico 
#      10.00000       10.00000       10.00000       25.00000       25.05044 
#  Leonardo.50%            GKi 
#      25.00000       25.00000

在这两种情况下,与median(rep(x, w))返回的结果相比,只有 isotoneGKi 给予相同的结果。

iaqfqrcu

iaqfqrcu8#

如果您正在使用survey包,假设您已定义调查设计,并且x是您感兴趣的变量:

svyquantile(~x, mydesign, c(0.5))
falq053o

falq053o9#

我来到这里是为了寻找加权分位数,所以我想我最好把我最终得到的留给未来的读者。自然,使用probs = 0.5将返回加权中位数。
我从MichaelChirico的answer开始,不幸的是它在边缘关闭,然后我决定从density()切换到approx(),最后,我相信我确定了校正因子,以确保与未加权quantile()的默认算法一致。

weighted.quantile <- function(x, w, probs = seq(0, 1, 0.25),
                              na.rm = FALSE, names = TRUE) {

  if (any(probs > 1) | any(probs < 0)) stop("'probs' outside [0,1]")

  if (length(w) == 1) w <- rep(w, length(x))
  if (length(w) != length(x)) stop("w must have length 1 or be as long as x")

  if (isTRUE(na.rm)) {
    w <- x[!is.na(x)]
    x <- x[!is.na(x)]
  }

  w <- w[order(x)] / sum(w)
  x <- x[order(x)]

  cum_w <- cumsum(w) - w * (1 - (seq_along(w) - 1) / (length(w) - 1))
  res <- approx(x = cum_w, y = x, xout = probs)$y

  if (isTRUE(names)) {
    res <- setNames(res, paste0(format(100 * probs, digits = 7), "%"))
  }
  res
}

权重统一时,加权分位数与常规未加权分位数相同:

x <- rnorm(100)
stopifnot(stopifnot(identical(weighted.quantile(x, w = 1), quantile(x)))

示例使用与weighted.mean()手册页中相同的数据。

x <- c(3.7, 3.3, 3.5, 2.8)
w <- c(5,   5,   4,   1)/15
stopifnot(isTRUE(all.equal(
  weighted.quantile(x, w, 0:4/4, names = FALSE),
  c(2.8, 3.33611111111111, 3.46111111111111, 3.58157894736842,
    3.7)
)))

这是给那些只想得到加权中值的人的:

weighted.median <- function(x, w, ...) {
  weighted.quantile(x, w, probs = 0.5, names = FALSE, ...)
}

相关问题