如何找到统计模式？

Question

R中，mean() 和 median() 是标准函数，它们执行你期望的操作。 mode() 返回对象内部存储模式，而不是其参数中出现最多的值。但是否有一个标准库函数可以实现向量（或列表）的统计众数？

- Nick

---

请问您的数据是整数、数值型还是分类型的？数值型的众数估计方法与其他类型不同，需要使用区间。请参考modeest。 - smci

为什么 R 语言没有内置的众数函数？为什么 R 认为 mode 和函数 class 是相同的？ - Corey Levinson

Answer 1

还有一种解决方法，适用于数值和字符/因子数据：

Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

在我的小型计算机上，可以在约半秒钟内生成并找到一个包含1000万个整数的向量的众数。

如果您的数据集可能有多个众数，则上述解决方案采用与which.max相同的方法，并返回众数集合中出现的第一个值。 要返回所有众数，请使用此变体（来自评论中的@digEmAll）：

Modes <- function(x) {
  ux <- unique(x)
  tab <- tabulate(match(x, ux))
  ux[tab == max(tab)]
}

- Ken Williams

---

也适用于逻辑向量！对于所有类型的向量都保留数据类型（不像其他答案中的某些实现）。 - DavidC

在多模态数据集的情况下（例如 c(1,1,2,2)），这并不返回所有模式。你应该将最后一行改为：tab <- tabulate(match(x, ux)); ux[tab == max(tab)]。 - digEmAll

为此，您需要将ux[which.max(tabulate(match(x, ux)))]替换为max(tabulate(match(x, ux)))。 - Ken Williams

注意到 Mode(1:3) 返回 1，而 Mode(3:1) 返回 3，因此 Mode 函数返回最常见的元素，如果所有元素都是唯一的，则返回第一个元素。 - Enrique Pérez Herrero

正如恩里克所说：当没有众数时，这种方法会失败，而且会给你留下第一个值就是众数的印象。如果在这些情况下返回0或NA将会更好。 - not2qubit

Answer 2

在 R 邮件列表中找到这个，希望它有所帮助。这也是我一直在想的。您需要对数据进行 table()，排序然后选择第一个名称。这很 hackish 但应该能用。

names(sort(-table(x)))[1]

- Dan

---

这是一个巧妙的解决方法。它有一些缺点：排序算法可能比基于max()的方法更消耗空间和时间（因此在较大的样本列表中要避免使用）。此外，输出的模式是 "字符" 而不是 "数值"。（请原谅双关语/歧义）当然，需要测试多模分布通常需要存储排序表以避免重新计算。 - mjv

我使用了一个包含1百万个元素的因子进行运行时间测量，这个解决方案比被接受的答案快了近3倍！ - vonjd

我刚刚使用as.numeric()将它转换为数字。完美地运行了。谢谢！ - Abhishek Singh

这个解决方案的问题在于，在存在多个众数的情况下，它是不正确的。 - vonjd

Answer 3

有一个叫做modeest的包，它提供了单一峰值（有时也是多峰值）数据的估计模式以及通常概率分布的模式值。

mySamples <- c(19, 4, 5, 7, 29, 19, 29, 13, 25, 19)

library(modeest)
mlv(mySamples, method = "mfv")

Mode (most likely value): 19 
Bickel's modal skewness: -0.1 
Call: mlv.default(x = mySamples, method = "mfv")

更多信息请参见此页面

您还可以在CRAN任务视图：概率分布中寻找“模式估计”。已提出两个新软件包。

- George Dontas

---

要仅获取众数值，可以使用mfv(mySamples)[1]。这里的数字1很重要，因为它实际上返回了最频繁出现的值们。 - atomicules

这个例子好像不起作用： library(modeest) a <- rnorm(50, 30, 2) b <- rnorm(100, 35, 2) c <- rnorm(20, 37, 2) temperatureºC <- c(a, b, c) hist(temperatureºC) #平均值 abline(v=mean(temperatureºC),col="red",lwd=2) #中位数 abline(v=median(temperatureºC),col="black",lwd=2) #众数 abline(v=mlv(temperatureºC, method = "mfv")[1],col="orange",lwd=2) - Agus camacho

@atomicules：使用[1]只能得到第一个模式。对于双峰或一般的n模态分布，您只需要使用mfv（mySamples）。 - petzi

针对R 3.6.0版本，当我使用函数'mlv'时出现“could not find function 'mlv'”的错误信息，并且在尝试使用mfv(mysamples)时同样出现相同的错误信息。它是否已经过时了？ - Dr Nisha Arora

@DrNishaArora：你下载了'modeest'包吗？ - petzi

Answer 4

我认为Ken Williams在上面的帖子中写得很好，我添加了几行以处理NA值并将它变成了一个函数以便于使用。

我发现Ken Williams在上面的帖子中写得很好，我添加了几行代码来解决缺失值问题，并将其转换为一个函数，以方便使用。

Mode <- function(x, na.rm = FALSE) {
  if(na.rm){
    x = x[!is.na(x)]
  }

  ux <- unique(x)
  return(ux[which.max(tabulate(match(x, ux)))])
}

- jprockbelly

---

我已经找到了一些加速方法，请见下面的答案。 - Dan Houghton

Answer 5

一种快速而简单的估算连续单变量分布（例如正态分布）中的向量模式的方法是定义并使用以下函数：

estimate_mode <- function(x) {
  d <- density(x)
  d$x[which.max(d$y)]
}

然后获取模式估计：

x <- c(5.8, 5.6, 6.2, 4.1, 4.9, 2.4, 3.9, 1.8, 5.7, 3.2)
estimate_mode(x)
## 5.439788

- Rasmus Bååth

---

关于这个问题，只需要注意一点：你可以通过这种方法获得任何一组连续数字的“众数”（mode）。数据不需要来自正态分布也可以使用。以下是一个例子，使用均匀分布中的数字。set.seed(1);a<-runif(100);mode<-density(a)$x[which.max(density(a)$y)];abline(v=mode) - Jota

error in density.default(x, from = from, to = to) : need at least 2 points to select a bandwidth automatically - Sergio

@xhie 那个错误信息告诉你所有你需要知道的。如果你只有一个数据点，那么在调用 density 时需要手动设置带宽。然而，如果你只有一个数据点，那么该数据点的值很可能就是你对于众数的最佳猜测... - Rasmus Bååth

你说得没错，但我只加了一个小调整：estimate_mode <- function(x) { if (length(x)>1){ d <- density(x) d$x[which.max(d$y)] }else{ x } } 我正在测试一种估计主导风向的方法，而不是使用圆形包中的矢量平均方向。我正在处理多边形等级上的点，因此有时只有一个具有方向的点。谢谢！ - Sergio

@xhie 听起来很合理 :) - Rasmus Bååth

Answer 6

以下函数有三种形式：

方法 = "mode" [默认值]：计算单峰向量的众数，否则返回NA
方法 = "nmodes"：计算向量中的模数数量
方法 = "modes"：列出单峰或多模向量的所有模数。

modeav <- function (x, method = "mode", na.rm = FALSE)
{
  x <- unlist(x)
  if (na.rm)
    x <- x[!is.na(x)]
  u <- unique(x)
  n <- length(u)
  #get frequencies of each of the unique values in the vector
  frequencies <- rep(0, n)
  for (i in seq_len(n)) {
    if (is.na(u[i])) {
      frequencies[i] <- sum(is.na(x))
    }
    else {
      frequencies[i] <- sum(x == u[i], na.rm = TRUE)
    }
  }
  #mode if a unimodal vector, else NA
  if (method == "mode" | is.na(method) | method == "")
  {return(ifelse(length(frequencies[frequencies==max(frequencies)])>1,NA,u[which.max(frequencies)]))}
  #number of modes
  if(method == "nmode" | method == "nmodes")
  {return(length(frequencies[frequencies==max(frequencies)]))}
  #list of all modes
  if (method == "modes" | method == "modevalues")
  {return(u[which(frequencies==max(frequencies), arr.ind = FALSE, useNames = FALSE)])}  
  #error trap the method
  warning("Warning: method not recognised.  Valid methods are 'mode' [default], 'nmodes' and 'modes'")
  return()
}

- Chris

---

在您对这些函数的描述中，您交换了“modes”和“nmodes”。请查看代码。实际上，“nmodes”返回值向量，“modes”返回模式数。尽管如此，您的函数是我迄今为止见过的找到模式的最佳解决方案。 - Grzegorz Adam Kowalski

非常感谢您的评论。现在，“nmode”和“modes”应该按预期运行。 - Chris

你的函数几乎可以工作，除了当每个值使用 method = 'modes' 出现相等的次数时。然后函数返回所有唯一的值，但实际上没有众数，所以应该返回 NA。我会添加另一个答案，其中包含稍微优化过的函数版本，感谢你的启发！ - hugovdberg

这个函数通常只有在对多峰向量使用默认方法时，才会生成非空数值向量的NA。例如，像1、2、3、4这样的简单数字序列的模式实际上是序列中的所有数字，因此对于类似的序列，“modes”表现得符合预期。 例如，modeave(c(1,2,3,4), method = "modes") 返回 [1] 1 2 3 4 尽管如此，我仍然非常希望看到该函数得到优化，因为它在当前状态下相当耗费资源。 - Chris

对于这个函数的更高效版本，请参见上面@hugovdberg的帖子 :) - Chris

Answer 7

这里，另一个解决方案：

freq <- tapply(mySamples,mySamples,length)
#or freq <- table(mySamples)
as.numeric(names(freq)[which.max(freq)])

- teucer

---

你可以将第一行替换为table。 - Jonathan Chang

我认为“tapply”比“table”更有效率，但它们都使用for循环。我认为使用“table”的解决方案是等效的。我更新了答案。 - teucer

Answer 8

一份对 Ken Williams 答案的小修改，加上可选参数 na.rm 和 return_multiple。

与依赖于 names() 的答案不同，此答案在返回值中保持了 x 的数据类型。

stat_mode <- function(x, return_multiple = TRUE, na.rm = FALSE) {
  if(na.rm){
    x <- na.omit(x)
  }
  ux <- unique(x)
  freq <- tabulate(match(x, ux))
  mode_loc <- if(return_multiple) which(freq==max(freq)) else which.max(freq)
  return(ux[mode_loc])
}

为了展示它如何与可选参数配合并保持数据类型：

foo <- c(2L, 2L, 3L, 4L, 4L, 5L, NA, NA)
bar <- c('mouse','mouse','dog','cat','cat','bird',NA,NA)

str(stat_mode(foo)) # int [1:3] 2 4 NA
str(stat_mode(bar)) # chr [1:3] "mouse" "cat" NA
str(stat_mode(bar, na.rm=T)) # chr [1:2] "mouse" "cat"
str(stat_mode(bar, return_mult=F, na.rm=T)) # chr "mouse"

感谢@Frank进行简化处理。

- C8H10N4O2

Answer 9

基于 @Chris 的函数来计算众数或相关指标，然而使用 Ken Williams 的方法来计算频率。这个函数为没有众数的情况提供了修复（所有元素同样频繁），并且提供了一些更易读的 方法 名称。

Mode <- function(x, method = "one", na.rm = FALSE) {
  x <- unlist(x)
  if (na.rm) {
    x <- x[!is.na(x)]
  }

  # Get unique values
  ux <- unique(x)
  n <- length(ux)

  # Get frequencies of all unique values
  frequencies <- tabulate(match(x, ux))
  modes <- frequencies == max(frequencies)

  # Determine number of modes
  nmodes <- sum(modes)
  nmodes <- ifelse(nmodes==n, 0L, nmodes)

  if (method %in% c("one", "mode", "") | is.na(method)) {
    # Return NA if not exactly one mode, else return the mode
    if (nmodes != 1) {
      return(NA)
    } else {
      return(ux[which(modes)])
    }
  } else if (method %in% c("n", "nmodes")) {
    # Return the number of modes
    return(nmodes)
  } else if (method %in% c("all", "modes")) {
    # Return NA if no modes exist, else return all modes
    if (nmodes > 0) {
      return(ux[which(modes)])
    } else {
      return(NA)
    }
  }
  warning("Warning: method not recognised.  Valid methods are 'one'/'mode' [default], 'n'/'nmodes' and 'all'/'modes'")
}

由于它使用Ken的方法来计算频率，因此性能也得到了优化。使用AkselA的帖子作为基准，我对之前的一些答案进行了基准测试，以展示我的函数在性能上与Ken的函数接近，各种输出选项的条件语句仅会导致轻微的开销：

- hugovdberg

---

我刚刚查看了一下，但是你指的是哪个版本的 pracma 包？据我所知，1.9.3 版本有完全不同的实现。 - hugovdberg

函数的修改很好。经过进一步阅读，我得出结论：关于均匀分布或单频分布是否具有节点，没有共识，一些来源说模式列表是分布本身，而其他人则说没有节点。唯一的共识是为这样的分布生成模式列表既不是非常信息丰富，也不是特别有意义的。 如果您希望该函数在这种情况下产生模式，则删除以下行：nmodes <- ifelse(nmodes==n, 0L, nmodes) - Chris

@ greendiod 抱歉，我错过了你的评论。您可以通过此Gist获取：https://gist.github.com/Hugovdberg/0f00444d46efd99ed27bbe227bdc4d37 - hugovdberg

这可能是最强大的答案！ - not2qubit

我建议在函数结尾处使用 stop 而不是 warning，因为该函数没有任何有意义的返回值。 - Gregor Thomas

Answer 10

collapse包中的通用函数fmode现已在CRAN上提供，它实现了基于索引哈希的C++模式。它比上述任何方法都要快得多。它附带有向量、矩阵、数据框和dplyr分组tibbles的方法。语法：

library(collapse)
fmode(x, g = NULL, w = NULL, ...)

其中x可以是上述对象之一，g提供一个可选的分组向量或分组向量列表（用于在C++中执行分组模式计算），w（可选）提供一个数字权重向量。在分组tibble方法中，没有g参数，您可以使用data %>% group_by(idvar) %>% fmode。

- Sebastian