计算 R 中多峰分布中的模式答案

【问题标题】：Calculate the modes in a multimodal distribution in R计算 R 中多峰分布中的模式
【发布时间】：2015-02-09 16:15:59
【问题描述】：

我测量了我所有孩子的身高。当我沿长度轴绘制所有高度时，结果如下所示：

每个红色（男孩）或紫色（女孩）勾号都是一个孩子。如果两个孩子的身高相同（以毫米为单位），蜱就会堆积起来。目前有七个身高相同的孩子。（刻度的高度和宽度没有意义。它们已被缩放为可见。）

如您所见，不同的高度不是沿轴均匀分布，而是围绕某些值聚集。

数据的直方图和密度图如下所示（在this answer 之后绘制了两个密度估计值）：

如您所见，这是一个多模态分布。

如何计算模式（R 中）？

以下是供您使用的原始数据：

mm <- c(418, 527, 540, 553, 554, 558, 613, 630, 634, 636, 645, 648, 708, 714, 715, 725, 806, 807, 822, 823, 836, 837, 855, 903, 908, 910, 911, 913, 915, 923, 935, 945, 955, 957, 958, 1003, 1006, 1015, 1021, 1021, 1022, 1034, 1043, 1048, 1051, 1054, 1058, 1100, 1102, 1103, 1117, 1125, 1134, 1138, 1145, 1146, 1150, 1152, 1210, 1211, 1213, 1223, 1226, 1334)

【问题讨论】：

我想我很困惑。你想找到mm的模式吗？它不是多模态向量，因为众数是 1021。您需要使用 mm 进行任何事先计算吗？
@LyzadeR 请参阅en.wikipedia.org/wiki/Multimodal_distribution 简化版：我想要的是我的问题中密度曲线的顶点。
孩子太多了。你的后宫有多大？
@ChuckD。有趣的。有人对此发表评论花了 5 个小时 :-)
从单峰分布中提取的少量观察结果通常如下所示：（例如，set.seed(6); hist(rnorm(64))），并且“找到模式”的过程定义不明确，因为有很多方法这是可以做到的……

标签： r distribution

【解决方案1】：

我使用你的 mm 数据自己构建了一些东西。

首先让我们绘制 mm 的密度以便可视化模式：

plot(density(mm))

因此，我们可以看到此分布中有 2 种模式。一个在 600 左右，一个在 1000 左右。让我们看看如何找到它们。

为了找到模式索引，我做了这个函数：

find_modes<- function(x) {
  modes <- NULL
  for ( i in 2:(length(x)-1) ){
    if ( (x[i] > x[i-1]) & (x[i] > x[i+1]) ) {
      modes <- c(modes,i)
    }
  }
  if ( length(modes) == 0 ) {
    modes = 'This is a monotonic distribution'
  }
  return(modes)
}

让我们试试我们的密度：

mymodes_indices <- find_modes(density(mm)$y) #you need to try it on the y axis

现在mymodes_indices 包含我们模式的索引，即：

> density(mm)$y[mymodes_indices]  #just to confirm that those are the correct
[1] 0.0008946929 0.0017766183

> density(mm)$x[mymodes_indices] #the actual modes
[1]  660.2941 1024.9067

希望对你有帮助！

【讨论】：

此解决方案高度依赖于平滑参数和底层分布，不是解决此问题的通用方法，并且会在噪声set.seed(6); find_modes(density(runif(64))$y) 中找到“模式”

【解决方案2】：

我修改了Peak of the kernel density estimation 中Jeffrey Evans 的答案，以允许修改bw 参数并相应地获得或多或少的峰值。对于其他情况，这是必要的，这将产生许多具有公认答案的峰值。参数signifi 允许处理关系。

library(dplyr)
library(tidyr)
get.modes2 <- function(x,adjust,signifi,from,to) {  
  den <- density(x, kernel=c("gaussian"),adjust=adjust,from=from,to=to)
  den.s <- smooth.spline(den$x, den$y, all.knots=TRUE, spar=0.1)
  s.1 <- predict(den.s, den.s$x, deriv=1)
  s.0 <- predict(den.s, den.s$x, deriv=0)
  den.sign <- sign(s.1$y)
  a<-c(1,1+which(diff(den.sign)!=0))
  b<-rle(den.sign)$values
  df<-data.frame(a,b)
  df = df[which(df$b %in% -1),]
  modes<-s.1$x[df$a]
  density<-s.0$y[df$a]
  df2<-data.frame(modes,density)
  df2$sig<-signif(df2$density,signifi)
  df2<-df2[with(df2, order(-sig)), ] 
  #print(df2)
  df<-as.data.frame(df2 %>% 
                      mutate(m = min_rank(desc(sig)) ) %>% #, count = sum(n)) %>% 
                      group_by(m) %>% 
                      summarize(a = paste(format(round(modes,2),nsmall=2), collapse = ',')) %>%
                      spread(m, a, sep = ''))
  colnames(df)<-paste0("m",1:length(colnames(df)))
  print(df)
}
mm <- c(418, 527, 540, 553, 554, 558, 613, 630, 634, 636, 645, 648, 708, 714, 715, 725, 806, 807, 822, 823, 836, 837, 855, 903, 908, 910, 911, 913, 915, 923, 935, 945, 955, 957, 958, 1003, 1006, 1015, 1021, 1021, 1022, 1034, 1043, 1048, 1051, 1054, 1058, 1100, 1102, 1103, 1117, 1125, 1134, 1138, 1145, 1146, 1150, 1152, 1210, 1211, 1213, 1223, 1226, 1334)
mmdf<-data.frame(mm=mm)
library(ggplot2)
#0.25 defines the number of peaks.
ggplot(mmdf,aes(mm)) + geom_density(adjust=0.25) + xlim((min(mm)-1),(max(mm)+1) )
#2 defines ties
modes<-get.modes2(mm,adjust=0.25,2,min(mm)-1,max(mm)+1)
#       m1     m2      m3            m4      m5     m6     m7              m8
#1 1031.40 921.81 1133.79 636.17,826.60 1216.43 548.14 715.22  418.80,1335.00

【讨论】：