如何找到统计模式?

在R中,mean()median()是你所期望的标准函数。mode()告诉你对象的内部存储模式,而不是参数中出现次数最多的值。但是是否存在一个标准库函数来实现向量(或列表)的统计模式?

366199 次浏览

在r邮件列表中发现了这个,希望对你有帮助。我也是这么想的。您将希望table()数据,排序,然后选择第一个名称。这有点粗俗,但应该有用。

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

R有如此多的附加包,其中一些可以很好地提供数字列表/系列/向量的[统计]模式。

然而,R的标准库本身似乎没有这样一个内置的方法!解决这个问题的一种方法是使用一些像下面这样的结构(如果你经常使用…则将其转换为函数):

mySamples <- c(19, 4, 5, 7, 29, 19, 29, 13, 25, 19)
tabSmpl<-tabulate(mySamples)
SmplMode<-which(tabSmpl== max(tabSmpl))
if(sum(tabSmpl == max(tabSmpl))>1) SmplMode<-NA
> SmplMode
[1] 19

对于更大的示例列表,应该考虑使用一个临时变量max(tabSmpl)值(我不知道R会自动优化这个)

参考:参见KickStarting R lesson .
中的“How about median and mode? 这似乎证实了(至少在写这节课的时候)R中没有模态函数(嗯…

. Mode(),正如你所发现的,用于断言变量的类型)

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任务视图:概率分布。已经提出了两个新的一揽子计划。

这里有另一个解决方案:

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

为了生成模式,我写了下面的代码。

MODE <- function(dataframe){
DF <- as.data.frame(dataframe)


MODE2 <- function(x){
if (is.numeric(x) == FALSE){
df <- as.data.frame(table(x))
df <- df[order(df$Freq), ]
m <- max(df$Freq)
MODE1 <- as.vector(as.character(subset(df, Freq == m)[, 1]))


if (sum(df$Freq)/length(df$Freq)==1){
warning("No Mode: Frequency of all values is 1", call. = FALSE)
}else{
return(MODE1)
}


}else{
df <- as.data.frame(table(x))
df <- df[order(df$Freq), ]
m <- max(df$Freq)
MODE1 <- as.vector(as.numeric(as.character(subset(df, Freq == m)[, 1])))


if (sum(df$Freq)/length(df$Freq)==1){
warning("No Mode: Frequency of all values is 1", call. = FALSE)
}else{
return(MODE1)
}
}
}


return(as.vector(lapply(DF, MODE2)))
}

让我们试试吧:

MODE(mtcars)
MODE(CO2)
MODE(ToothGrowth)
MODE(InsectSprays)

还有一个解决方案,它适用于数字&字符/因素数据:

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

在我的小机器上,它可以生成&在半秒内找到一个10m整型向量的模态。

如果你的数据集可能有多种模式,上面的解决方案采用与which.max相同的方法,并返回模式集的首次出现值。要返回所有模式,使用这个变体(来自评论中的@digEmAll):

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

另一个按频率排序的简单选项是使用rle:

df = as.data.frame(unclass(rle(sort(mySamples))))
df = df[order(-df$lengths),]
head(df)

估计来自连续单变量分布(例如正态分布)的数字向量的模式的一种快速而肮脏的方法是定义并使用以下函数:

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

下面的函数有三种形式:

method = "mode"[默认值]:计算单模态向量的模式,否则返回NA
Method = "nmodes":计算向量
中的模式数 Method = "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()
}

抱歉,我可能把它理解得太简单了,但这不是可以工作的吗?(我的机器上的1E6值在1.3秒内):

t0 <- Sys.time()
summary(as.factor(round(rnorm(1e6), 2)))[1]
Sys.time()-t0

你只需要用你的向量替换“round(rnorm(1e6),2)”。

我还不能投票,但Rasmus Bååth的答案是我正在寻找的。 但是,我将稍微修改一下,允许将分布限制在0到1之间。< / p >

estimate_mode <- function(x,from=min(x), to=max(x)) {
d <- density(x, from=from, to=to)
d$x[which.max(d$y)]
}

我们知道你可能根本不想约束你的分布,那么设置from=-"BIG NUMBER", to="BIG NUMBER"

您还可以计算一个实例在您的集合中出现的次数,并找到最大次数。如。

> temp <- table(as.vector(x))
> names (temp)[temp==max(temp)]
[1] "1"
> as.data.frame(table(x))
r5050 Freq
1     0   13
2     1   15
3     2    6
>

效果很好

> a<-c(1,1,2,2,3,3,4,4,5)
> names(table(a))[table(a)==max(table(a))]

可以尝试以下功能:

  1. 将数值转换为因子
  2. 使用summary()获取频率表
  3. 返回模式频率最大的索引
  4. 转换因子回到数字,即使有超过1个模式,这个函数工作良好!
mode <- function(x){
y <- as.factor(x)
freq <- summary(y)
mode <- names(freq)[freq[names(freq)] == max(freq)]
as.numeric(mode)
}

我将使用density()函数来确定一个(可能是连续的)分布的平滑最大值:

function(x) density(x, 2)$x[density(x, 2)$y == max(density(x, 2)$y)]

其中x是数据集合。注意密度函数的调整参数,它调节平滑。

我发现Ken Williams上面的帖子很棒,我添加了几行来解释NA值,并使其成为一个函数。

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)))])
}

虽然我喜欢肯威廉姆斯简单的功能,我想检索多种模式,如果他们存在。考虑到这一点,我使用下面的函数,它返回多个模式或单个模式的列表。

rmode <- function(x) {
x <- sort(x)
u <- unique(x)
y <- lapply(u, function(y) length(x[x==y]))
u[which( unlist(y) == max(unlist(y)) )]
}

下面是一个查找模式的函数:

mode <- function(x) {
unique_val <- unique(x)
counts <- vector()
for (i in 1:length(unique_val)) {
counts[i] <- length(which(x==unique_val[i]))
}
position <- c(which(counts==max(counts)))
if (mean(counts)==max(counts))
mode_x <- 'Mode does not exist'
else
mode_x <- unique_val[position]
return(mode_x)
}

另一个可能的解决方案:

Mode <- function(x) {
if (is.numeric(x)) {
x_table <- table(x)
return(as.numeric(names(x_table)[which.max(x_table)]))
}
}

用法:

set.seed(100)
v <- sample(x = 1:100, size = 1000000, replace = TRUE)
system.time(Mode(v))

输出:

   user  system elapsed
0.32    0.00    0.31

我浏览了所有这些选项,开始想知道它们的相对特性和性能,所以我做了一些测试。如果其他人也好奇,我在这里分享我的结果。

我不想为这里发布的所有函数而烦恼,我选择了一个基于一些标准的示例:函数应该对字符、因子、逻辑和数字向量都有效,它应该适当地处理na和其他有问题的值,输出应该是“合理的”,即没有数字作为字符或其他类似的愚蠢行为。

我还添加了一个我自己的函数,它基于与chrispy相同的rle思想,除了用于更一般的用途:

library(magrittr)


Aksel <- function(x, freq=FALSE) {
z <- 2
if (freq) z <- 1:2
run <- x %>% as.vector %>% sort %>% rle %>% unclass %>% data.frame
colnames(run) <- c("freq", "value")
run[which(run$freq==max(run$freq)), z] %>% as.vector
}


set.seed(2)


F <- sample(c("yes", "no", "maybe", NA), 10, replace=TRUE) %>% factor
Aksel(F)


# [1] maybe yes


C <- sample(c("Steve", "Jane", "Jonas", "Petra"), 20, replace=TRUE)
Aksel(C, freq=TRUE)


# freq value
#    7 Steve

我最终通过microbenchmark在两组测试数据上运行了五个函数。函数名指的是它们各自的作者:

enter image description here

Chris的函数在默认情况下被设置为method="modes"na.rm=TRUE,以使其更具可比性,但除此之外,函数的作者在这里使用。

单就速度而言,Kens版本轻松获胜,但它也是唯一一个只报告一种模式的版本,不管有多少种模式。通常情况下,在速度和多功能性之间需要权衡。在method="mode"中,Chris的版本将返回一个值,如果有一个模式,否则为NA。我觉得这招不错。 我还认为,有趣的是,一些函数会受到不断增加的惟一值的影响,而另一些函数则几乎不受影响。我还没有详细研究代码来弄清楚为什么会这样,除了消除逻辑/数字的原因

基于@Chris的函数来计算模态或相关指标,但是使用Ken Williams的方法来计算频率。这个函数修复了完全没有模式的情况(所有元素的频率相同),并提供了一些更易读的method名称。

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的,各种输出选项的条件只导致很小的开销:  Mode函数的比较 < / p >

计算包含离散值的向量“v”的MODE的一个简单方法是:

names(sort(table(v)))[length(sort(table(v)))]

这个黑客应该工作良好。给你的值以及模式的计数:

Mode <- function(x){
a = table(x) # x is a vector
return(a[which.max(a)])
}

计算模式大多是在有因素变量的情况下才可以使用

labels(table(HouseVotes84$V1)[as.numeric(labels(max(table(HouseVotes84$V1))))])

HouseVotes84是在“mlbench”包中可用的数据集。

它会给出最大标签值。它更容易由内置函数本身使用,而无需编写函数。

下面是可以用来找到R中矢量变量的模式的代码。

a <- table([vector])


names(a[a==max(a)])

对Ken Williams的回答做了一个小修改,增加了可选参数na.rmreturn_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的简化。

对此有多种解决方案。我检查了第一个,然后写了我自己的。把它贴在这里,如果它能帮助到任何人:

Mode <- function(x){
y <- data.frame(table(x))
y[y$Freq == max(y$Freq),1]
}

让我们用几个例子来测试一下。我正在使用iris数据集。让我们用数值数据进行测试

> Mode(iris$Sepal.Length)
[1] 5

你可以验证这是正确的。

现在虹膜数据集中唯一的非数字字段(Species)没有模式。让我们用我们自己的例子进行测试

> test <- c("red","red","green","blue","red")
> Mode(test)
[1] red

编辑

正如注释中提到的,用户可能希望保留输入类型。在这种情况下,mode函数可以修改为:

Mode <- function(x){
y <- data.frame(table(x))
z <- y[y$Freq == max(y$Freq),1]
as(as.character(z),class(x))
}

函数的最后一行只是将最终的模式值强制为原始输入的类型。

模式并不是在所有情况下都有用。所以函数应该处理这种情况。试试下面的函数。

Mode <- function(v) {
# checking unique numbers in the input
uniqv <- unique(v)
# frquency of most occured value in the input data
m1 <- max(tabulate(match(v, uniqv)))
n <- length(tabulate(match(v, uniqv)))
# if all elements are same
same_val_check <- all(diff(v) == 0)
if(same_val_check == F){
# frquency of second most occured value in the input data
m2 <- sort(tabulate(match(v, uniqv)),partial=n-1)[n-1]
if (m1 != m2) {
# Returning the most repeated value
mode <- uniqv[which.max(tabulate(match(v, uniqv)))]
} else{
mode <- "Two or more values have same frequency. So mode can't be calculated."
}
} else {
# if all elements are same
mode <- unique(v)
}
return(mode)
}

输出,

x1 <- c(1,2,3,3,3,4,5)
Mode(x1)
# [1] 3


x2 <- c(1,2,3,4,5)
Mode(x2)
# [1] "Two or more varibles have same frequency. So mode can't be calculated."


x3 <- c(1,1,2,3,3,4,5)
Mode(x3)
# [1] "Two or more values have same frequency. So mode can't be calculated."

这建立在jprockbelly的答案上,通过对非常短的向量增加速度。这在将mode应用到data.frame或包含很多小组的数据表时非常有用:

Mode <- function(x) {
if ( length(x) <= 2 ) return(x[1])
if ( anyNA(x) ) x = x[!is.na(x)]
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}

我假设你的观测值是来自实数,当你的观测值是2,2,3,3时,你期望模式为2.5,然后你可以用mode = l1 + i * (f1-f0) / (2f1 - f0 - f2)估计模式,其中l1..最频繁类的下限,f1..最频繁类的频率,f0..在最频繁类之前的类的频率,f2..在最频繁类之后的类的频率和..给出的类间隔,例如在10, 1中:

#Small Example
x <- c(2,2,3,3) #Observations
i <- 1          #Class interval


z <- hist(x, breaks = seq(min(x)-1.5*i, max(x)+1.5*i, i), plot=F) #Calculate frequency of classes
mf <- which.max(z$counts)   #index of most frequent class
zc <- z$counts
z$breaks[mf] + i * (zc[mf] - zc[mf-1]) / (2*zc[mf] - zc[mf-1] - zc[mf+1])  #gives you the mode of 2.5




#Larger Example
set.seed(0)
i <- 5          #Class interval
x <- round(rnorm(100,mean=100,sd=10)/i)*i #Observations


z <- hist(x, breaks = seq(min(x)-1.5*i, max(x)+1.5*i, i), plot=F)
mf <- which.max(z$counts)
zc <- z$counts
z$breaks[mf] + i * (zc[mf] - zc[mf-1]) / (2*zc[mf] - zc[mf-1] - zc[mf+1])  #gives you the mode of 99.5

如果你想要最频繁级别并且你有不止一个最频繁的级别,你可以得到所有它们,例如:

x <- c(2,2,3,5,5)
names(which(max(table(x))==table(x)))
#"2" "5"

在我看来,如果一个集合有一个模式,那么它的元素就可以与自然数一一对应。因此,查找模式的问题简化为生成这样一个映射,查找映射值的模式,然后映射回集合中的一些项。(处理NA发生在映射阶段)。

我有一个histogram函数,它对类似的原理进行操作。(在这里给出的代码中使用的特殊函数和操作符应该在夏皮罗和/或neatOveRse中定义。在此复制夏皮罗和奈尔斯的部分是经过允许的;复制的片段可根据本网站的条款使用。)R 伪代码 for histogram

.histogram <- function (i)
if (i %|% is.empty) integer() else
vapply2(i %|% max %|% seqN, `==` %<=% i %O% sum)


histogram <- function(i) i %|% rmna %|% .histogram

(特殊的二进制运算符完成管道局部套用作文)我还有一个maxloc函数,它类似于which.max,但返回所有一个向量的绝对最大值。R 伪代码 for maxloc

FUNloc <- function (FUN, x, na.rm=F)
which(x == list(identity, rmna)[[na.rm %|% index.b]](x) %|% FUN)


maxloc <- FUNloc %<=% max


minloc <- FUNloc %<=% min # I'M THROWING IN minloc TO EXPLAIN WHY I MADE FUNloc

然后

imode <- histogram %O% maxloc

而且

x %|% map %|% imode %|% unmap

将计算任何集合的模式,只要定义了适当的map-ping和unmap-ping函数。

添加raster::modal()作为一个选项,尽管注意raster是一个很大的包,如果你不做地理空间方面的工作,可能不值得安装。

对于那些特别热衷的人来说,源代码可以从https://github.com/rspatial/raster/blob/master/src/modal.cpphttps://github.com/rspatial/raster/blob/master/R/modal.R中提取到个人R包中。

现在CRAN上可用的collapse包中的泛型函数fmode实现了基于索引哈希的基于c++的模式。它比上述任何一种方法都要快得多。它提供了向量、矩阵、data.frames和dplyr分组tibbles的方法。语法:

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

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

如果你问R中的内置函数,也许你可以在包pracma中找到它。在这个包中,有一个名为Mode的函数。

这是我的数据。返回完整表的逐行模式的表解决方案。我用它来推断行类。它负责data中新的set()函数。桌子,应该很快。虽然它不管理NA,但可以通过查看本页上的众多其他解决方案添加。

majorityVote <- function(mat_classes) {
#mat_classes = dt.pour.centroids_num
dt.modes <- data.table(mode = integer(nrow(mat_classes)))
for (i in 1:nrow(mat_classes)) {
cur.row <- mat_classes[i]
cur.mode <- which.max(table(t(cur.row)))
set(dt.modes, i=i, j="mode", value = cur.mode)
}


return(dt.modes)
}

可能的用法:

newClass <- majorityVote(my.dt)  # just a new vector with all the modes