如何在分组数据集上计算中位数?

9
我的数据集如下:
salary  number
1500-1600   110
1600-1700   180
1700-1800   320
1800-1900   460
1900-2000   850
2000-2100   250
2100-2200   130
2200-2300   70
2300-2400   20
2400-2500   10

我该如何计算这个数据集的中位数?以下是我尝试过的方法:

x <- c(110, 180, 320, 460, 850, 250, 130, 70, 20, 10)
colnames <- "numbers"
rownames <- c("[1500-1600]", "(1600-1700]", "(1700-1800]", "(1800-1900]", 
              "(1900-2000]", "(2000,2100]", "(2100-2200]", "(2200-2300]",
              "(2300-2400]", "(2400-2500]")
y <- matrix(x, nrow=length(x), dimnames=list(rownames, colnames))
data.frame(y, "cumsum"=cumsum(y))

            numbers cumsum
[1500-1600]     110    110
(1600-1700]     180    290
(1700-1800]     320    610
(1800-1900]     460   1070
(1900-2000]     850   1920
(2000,2100]     250   2170
(2100-2200]     130   2300
(2200-2300]      70   2370
(2300-2400]      20   2390
(2400-2500]      10   2400

在这里,你可以看到中位数的频率为2400/2=1200。它在10701920之间。因此中位数类别(1900-2000]组。您可以使用下面的公式得出这个结果:

中位数 = L + h/f (n/2 - c)

其中:

L 是中位数类别的下限边界
h 是中位数类别的大小,即中位数类别的上下限边界之差
f 是中位数类别的频率
c 是中位数类别的前一个累计频率
n/2 是观察次数总数除以2(即总和f/2)

或者,通过以下方法定义中位数类别

在累计频率列中找到n/2。

获取该类别的区间。

代码如下:

> 1900 + (1200 - 1070) / (1920 - 1070) * (2000 - 1900)    
[1] 1915.294

现在我想使上述表达更加优雅 - 即1900+(1200-1070)/(1920-1070)*(2000-1900)。我该如何实现?

3
你尝试过什么?你考虑过提供一个可重现的例子吗?https://dev59.com/eG025IYBdhLWcg3whGSx - Roman Luštrik
6个回答

7

既然你已经知道了公式,那么创建一个函数来进行计算应该很容易。

这里我创建了一个基本的函数来帮助你入门。该函数需要四个参数:

  • frequencies:频率向量(在你的第一个例子中是“number”)
  • intervals:一个 2 行矩阵,与频率长度相同的列数,第一行是下限,第二行是上限。或者,“intervals”可以是你的 data.frame 中的一列,并且你可以指定 sep(可能还有 trim)来让函数自动为你创建所需的矩阵。
  • sep:你的 data.frame 中“intervals”列中的分隔符字符。
  • trim:需要在尝试强制转换为数字矩阵之前删除的字符的正则表达式。一个模式内置于函数中:trim = "cut"。这将设置正则表达式模式以从输入中删除 (、)、[ 和 ]。

以下是该函数(其中包含注释,显示我如何使用你的说明来组合它):

GroupedMedian <- function(frequencies, intervals, sep = NULL, trim = NULL) {
  # If "sep" is specified, the function will try to create the 
  #   required "intervals" matrix. "trim" removes any unwanted 
  #   characters before attempting to convert the ranges to numeric.
  if (!is.null(sep)) {
    if (is.null(trim)) pattern <- ""
    else if (trim == "cut") pattern <- "\\[|\\]|\\(|\\)"
    else pattern <- trim
    intervals <- sapply(strsplit(gsub(pattern, "", intervals), sep), as.numeric)
  }

  Midpoints <- rowMeans(intervals)
  cf <- cumsum(frequencies)
  Midrow <- findInterval(max(cf)/2, cf) + 1
  L <- intervals[1, Midrow]      # lower class boundary of median class
  h <- diff(intervals[, Midrow]) # size of median class
  f <- frequencies[Midrow]       # frequency of median class
  cf2 <- cf[Midrow - 1]          # cumulative frequency class before median class
  n_2 <- max(cf)/2               # total observations divided by 2

  unname(L + (n_2 - cf2)/f * h)
}

这是一个可供操作的样例 data.frame

mydf <- structure(list(salary = c("1500-1600", "1600-1700", "1700-1800", 
    "1800-1900", "1900-2000", "2000-2100", "2100-2200", "2200-2300", 
    "2300-2400", "2400-2500"), number = c(110L, 180L, 320L, 460L, 
    850L, 250L, 130L, 70L, 20L, 10L)), .Names = c("salary", "number"), 
    class = "data.frame", row.names = c(NA, -10L))
mydf
#       salary number
# 1  1500-1600    110
# 2  1600-1700    180
# 3  1700-1800    320
# 4  1800-1900    460
# 5  1900-2000    850
# 6  2000-2100    250
# 7  2100-2200    130
# 8  2200-2300     70
# 9  2300-2400     20
# 10 2400-2500     10

现在,我们可以简单地执行以下操作:
GroupedMedian(mydf$number, mydf$salary, sep = "-")
# [1] 1915.294

以下是一些虚构数据的函数示例:
set.seed(1)
x <- sample(100, 100, replace = TRUE)
y <- data.frame(table(cut(x, 10)))
y
#           Var1 Freq
# 1   (1.9,11.7]    8
# 2  (11.7,21.5]    8
# 3  (21.5,31.4]    8
# 4  (31.4,41.2]   15
# 5    (41.2,51]   13
# 6    (51,60.8]    5
# 7  (60.8,70.6]   11
# 8  (70.6,80.5]   15
# 9  (80.5,90.3]   11
# 10  (90.3,100]    6

### Here's GroupedMedian's output on the grouped data.frame...
GroupedMedian(y$Freq, y$Var1, sep = ",", trim = "cut")
# [1] 49.49231

### ... and the output of median on the original vector
median(x)
# [1] 49.5

顺便提一句,根据你提供的样本数据,我认为你的某个范围存在错误(所有范围都用短横线分隔,除了一个用逗号分隔),因为strsplit默认使用正则表达式进行分割,所以你可以这样使用该函数:
x<-c(110,180,320,460,850,250,130,70,20,10)
colnames<-c("numbers")
rownames<-c("[1500-1600]","(1600-1700]","(1700-1800]","(1800-1900]",
            "(1900-2000]"," (2000,2100]","(2100-2200]","(2200-2300]",
            "(2300-2400]","(2400-2500]")
y<-matrix(x,nrow=length(x),dimnames=list(rownames,colnames))
GroupedMedian(y[, "numbers"], rownames(y), sep="-|,", trim="cut")
# [1] 1915.294

似乎你的输入数据需要按升序排列。如果它们不按顺序排列,程序就不会对其进行必要的重新排列。 - oatmilkyway

4

我这样写是为了清楚地解释它的工作原理。下面附有更简洁的版本。

library(data.table)

#constructing the dataset with the salary range split into low and high
salarydata <- data.table(
  salaries_low = 100*c(15:24),
  salaries_high = 100*c(16:25),
  numbers = c(110,180,320,460,850,250,130,70,20,10)
)

#calculating cumulative number of observations
salarydata <- salarydata[,cumnumbers := cumsum(numbers)]
salarydata
   # salaries_low salaries_high numbers cumnumbers
   # 1:         1500          1600     110        110
   # 2:         1600          1700     180        290
   # 3:         1700          1800     320        610
   # 4:         1800          1900     460       1070
   # 5:         1900          2000     850       1920
   # 6:         2000          2100     250       2170
   # 7:         2100          2200     130       2300
   # 8:         2200          2300      70       2370
   # 9:         2300          2400      20       2390
   # 10:         2400          2500      10       2400

#identifying median group
mediangroup <- salarydata[
  (cumnumbers - numbers) <= (max(cumnumbers)/2) & 
  cumnumbers >= (max(cumnumbers)/2)]
mediangroup
   # salaries_low salaries_high numbers cumnumbers
   # 1:         1900          2000     850       1920

#creating the variables needed to calculate median
mediangroup[,l := salaries_low]
mediangroup[,h := salaries_high - salaries_low]
mediangroup[,f := numbers]
mediangroup[,c := cumnumbers- numbers]
n = salarydata[,sum(numbers)]

#calculating median
median <- mediangroup[,l + ((h/f)*((n/2)-c))]
median
   # [1] 1915.294

紧凑版 -

编辑:根据@AnandaMahto的建议改为函数。此外,使用更通用的变量名称。

library(data.table)

#Creating function

CalculateMedian <- function(
   LowerBound,
   UpperBound,
   Obs
)
{
   #calculating cumulative number of observations and n
   dataset <- data.table(UpperBound, LowerBound, Obs)

   dataset <- dataset[,cumObs := cumsum(Obs)]
   n = dataset[,max(cumObs)]

   #identifying mediangroup and dynamically calculating l,h,f,c. We already have n.
   median <- dataset[
      (cumObs - Obs) <= (max(cumObs)/2) & 
      cumObs >= (max(cumObs)/2),

      LowerBound + ((UpperBound - LowerBound)/Obs) * ((n/2) - (cumObs- Obs))
   ]

   return(median)
}


# Using function
CalculateMedian(
  LowerBound = 100*c(15:24),
  UpperBound = 100*c(16:25),
  Obs = c(110,180,320,460,850,250,130,70,20,10)
)
# [1] 1915.294

2
我个人希望看到你将你的答案转换为一个函数,就像我在我的答案中所做的那样。否则,这对于OP已经知道如何根据输入数据集手动计算答案来说并没有太大的增值。 - A5C1D2H2I1M1N2O1R2T1

3
(Sal <- sapply( strsplit(as.character(dat[[1]]), "-"), 
                                 function(x) mean( as.numeric(x) ) ) )
 [1] 1550 1650 1750 1850 1950 2050 2150 2250 2350 2450
require(Hmisc)
wtd.mean(Sal, weights = dat[[2]])
[1] 1898.75
wtd.quantile(Sal, weights=dat[[2]], probs=0.5)

一般化到加权中位数可能需要寻找具有此功能的软件包。

我猜你是指 weighted.mean 函数? - A5C1D2H2I1M1N2O1R2T1
1
不,我的意思是加权中位数。有些软件包确实有这样的功能。 - IRTFM
我的意思是在你的示例代码中。我不知道是否存在一个名为"wtd.mean"的现有函数 :) - A5C1D2H2I1M1N2O1R2T1
哦,对了。我总是加载rms/Hmisc,并且它们都有wtd.meanwtd.quantile函数,可以为0.5的分位数提供加权中位数。 - IRTFM

0
这种方法怎么样?为每个薪资档位创建向量,假设每个带宽均匀分布。然后从这些向量中创建一个大向量,并取中位数。与您类似,但结果略有不同。我不是数学家,所以这种方法可能不正确。
dat <- matrix(c(seq(1500, 2400, 100), seq(1600, 2500, 100), c(110, 180, 320, 460, 850, 250, 130, 70, 20, 10)), ncol=3)
median(unlist(apply(dat, 1, function(x) { ((1:x[3])/x[3])*(x[2]-x[1])+x[1] })))

返回1915.353


0
你尝试过如果是矩阵或数据框,使用medianapply(yourobject,2,median)吗?

-3

我认为这个概念应该适用于你。

$salaries = array(
       array("1500","1600"),
       array("1600","1700"),
       array("1700","1800"),
       array("1800","1900"),
       array("1900","2000"),
       array("2000","2100"),
       array("2100","2200"),
       array("2200","2300"),
       array("2300","2400"),
       array("2400","2500"),
      );
 $numbers = array("110","180","320","460","850","250","130","70","20","10");
 $cumsum = array();
 $n = 0;
 $count = 0;
 foreach($numbers as $key=>$number){    
$cumsum[$key] = $number;    
$n += $number;
if($count > 0){
    $cumsum[$key] += $cumsum[$key-1];       
}
++$count;
 }

 $classIndex = 0;
 foreach($cumsum as $key=>$cum){
if($cum < ($n/2)){
 $classIndex = $key+1;
}
 }
 $classRange = $salaries[$classIndex];
 $L = $classRange[0];
 $h = (float) $classRange[1] - $classRange[0];
 $f = $numbers[$classIndex];
 $c = $numbers[$classIndex-1];

 $Median = $L + ($h/$f)*(($n/2)-$c);
 echo $Median;

1
这是什么语言?OP正在寻找R语言的解决方案。 - A5C1D2H2I1M1N2O1R2T1
这是PHP,我想对于程序员来说概念已经足够了。 - Mohammad Ismail Khan
1
真的,但是将代码从一种语言翻译成另一种语言并不总是容易的,而在一种语言中可能非常高效的东西,在另一种语言中可能会变得非常慢。 - A5C1D2H2I1M1N2O1R2T1

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接