很久以前，我在地质文献中学到了一种有效的技术。（我很抱歉没有记住来源。）它包括研究内核密度估计器 (KDE) 的模式，因为带宽是变化的。

发生的情况是，在带宽非常大的情况下，数据看起来像一个单一模式的大块。这个使用 60 的带宽，其模式接近 110：

随着带宽的缩小，KDE 勾勒出人眼看到的更近的东西并出现更多的模式。这个使用 10 的带宽并具有三种明显的模式，第四种模式刚刚开始显示在 60 附近：

当带宽缩小太多时，KDE 就太详细了。这个带宽为 1 的有 36 种模式：

您可以使用“模式跟踪”来探索这种行为。对于整个范围内的每个带宽（从没有细节到太详细），它都会绘制模式。我跟踪了每种模式的演变并相应地为它们着色。例如，第一张图中的单模对应于中心红线（形状几乎像问号）；第二幅图中的四种模式对应于上升到高度（带宽）为 10 的 4 条迹线；第三幅图中的 36 种模式对应于所有 36 条迹线：

对带宽使用对数刻度可能是个好主意，如此处所示。

浏览模式跟踪将指示要识别的模式数量。 我选了四个。为了确定它们的位置，我发现所有带宽中迹线最垂直的点小于所有四种模式首次出现的位置：在这些位置，即使带宽发生变化，位置也是稳定的。令人欣慰的是（但不是真正必要的）所有四个位置都使用可比较的带宽。（如果沿迹线出现多个稳定点，确实应该更加小心：我会选择最大带宽小于所有模式出现的带宽的那个。）

找到模式后，我们可以将它们绘制在原始直方图上：

那么选择极限模式就很简单了。模式跟踪将告诉您它们的位置对您识别的模式数量和您使用的带宽有多敏感。在这个例子中，它表明最高模式在分裂成多个模式之前会随着更小的带宽变得更大，但其他三种模式保持相对稳定（它们的迹线在低带宽下几乎保持垂直）。

你选择什么形状的内核并不重要。 原始论文建议使用高斯核，我在这里做了。使用高斯并不等同于任何假设峰甚至近似具有高斯形状。因为高斯是（无限）平滑的，所以 KDE 也是如此，这意味着您可以使用微积分技术对它进行分析。

为了完全清楚，这里是模式迹线的数学说明。设核函数的单位面积和唯一模式为，数据为 带宽为的数据的 KDE 为卷积$K$$0$$x_1, \ldots, x_n.$$h\ge 0$

对于每个令为分布函数 数据的“模式轨迹”是的并集，因为的范围在区间上，其中被选择得如此之大，以至于包含所有$h\ge 0,$$M(h)$$x\to f(x,h).$$M(h)$$h$$(0, A)$$A$$M(h)$$h\ge A.$

模态迹线具有额外的结构：它可以分解（不一定是唯一的）为在区间上定义这种分解是最大的，因为任何两个不同的此类函数可能共有的唯一点是在它们域的端点处。我已经使用颜色来指定这些部分功能。$h$

除了选择要使用的模式数量（这在很大程度上取决于您对分析数据的正确分辨率的概念）之外，一切都可以自动化。 这是R我用来生成样本数据、分析它们并制作数字的代码。其结果将包含在记录模式跟踪的数据帧和包含有关所选模式的信息X的数组中。modes

顺便说一句，如果您自己编写代码，请注意使用快速傅里叶变换 (FFT) 最有效地获得 KDE。最有效的方法是对数据进行一次转换，然后将其乘以一系列转换后的内核，将每个乘积反转以生成 KDE。为了确定要搜索的带宽范围，使最大的大约为数据范围的四分之一，最小的可能是该范围的 3% 或 1%。

#
# Generate random values from a mixture distribution.
#
rmix <- function(n, mu, sigma, p) {
  matrix(rnorm(length(mu)*n, mu, sigma), ncol=n)[
         cbind(sample.int(length(mu), n, replace=TRUE, prob=p), 1:n)]
}
mu <- c(25, 60, 130, 190) # Means
sigma <- c(8, 13, 15, 19) # SDs
p <- c(.18, .2, .24, .28) # Relative proportions (needn't sum to 1)
n <- 1e4                  # Sample size
x <- rmix(n, mu, sigma, p)
#
# Find the modes of a KDE.
# (Quick and dirty: it assumes no mode spans more than one x value.)
#
findmodes <- function(kde) {
  kde$x[which(c(kde$y[-1],NA) < kde$y & kde$y > c(NA,kde$y[-length(kde$y)]))]
}
#
# Compute the mode trace by varying the bandwidth within a factor of 10 of
# the default bandwidth.  Track the modes as the bandwidth is decreased from
# its largest to its smallest value.
# This calculation is fast, so we can afford a detailed search.
#
m <- mean(x)
id <- 1
bw <- density(x)$bw * 10^seq(1,-1, length.out=101) 
modes.lst <- lapply(bw, function(h) {
  m.new <- sort(findmodes(density(x, bw=h)))
  # -- Associate each previous mode with a nearest new mode.
  if (length(m.new)==1) delta <- Inf else delta <- min(diff(m.new))/2
  d <- outer(m.new, m, function(x,y) abs(x-y))
  i <- apply(d, 2, which.min)
  g <- rep(NA_integer_, length(m.new))
  g[i] <- id[1:ncol(d)]
  #-- Create new ids for new modes that appear.
  k <- is.na(g)
  g[k] <- (sum(!k)+1):length(g)
  id <<- g
  m <<- m.new
  data.frame(bw=h, Mode=m.new, id=g)
})
X <- do.call(rbind, args=modes.lst)
X$id <- factor(X$id)
#
# Locate the modes at the most vertical portions of their traces.
#
minslope <- function(x, y) {
  f <- splinefun(x, y)
  e <- diff(range(x)) * 1e-4
  df2 <- function(x) ((f(x+e)-f(x-e)) / (2*e))^2 # Numerical derivative, squared
  v <- optimize(df2, c(min(x),max(x)))
  c(bw=v$minimum, slope=v$objective, Mode=f(v$minimum))
}
#
# Retain the desired modes.
#
n.modes <- 4 # USER SELECTED: Not automatic
bw.max <- max(subset(X, id==n.modes)$bw)
modes <- sapply(1:n.modes, function(i) {
  Y <- subset(X, id==i & bw <= bw.max)
  minslope(Y$bw, Y$Mode)
})
#
# Plot the results.
#
library(ggplot2)
ggplot(X, aes(bw, Mode)) +
  geom_line(aes(col=id), size=1.2, show.legend=FALSE) +
  geom_point(aes(bw, Mode), data=as.data.frame(t(modes)), size=3, col="Black", alpha=1/2) +
  scale_x_log10() +
  coord_flip() +
  ggtitle("Mode Trace")

ggplot(data.frame(x), aes(x, ..density..)) +
  geom_histogram(bins=500, fill="#2E75B2") +
  geom_vline(data=as.data.frame(t(modes)),
             mapping=aes(xintercept=Mode), col="#D18A4e", size=1) +
  ggtitle("Histogram With Modes")

这是一些多模态的假模拟数据。该图显示了三种图（用 R 制作）：（1）直方图 [蓝色]，（2）轴下方的刻度线，以及（3）数据的核密度估计器（KDE）[红色]。

hist(x, prob=T, br=40, col="skyblue2"); rug(x)
 lines(density(x), col="red")

我认为您的情节的主要困难在于它将本地联系与整体模式混淆了。我的直方图条的高度也不规则，因为我（故意）使用了太多条（通过使用参数br=40）。然而，即使选择了最佳的条数，直方图也不是寻找模式的最佳设备。

KDE 的“带宽”可以针对查找模式进行调整。上面，我使用了默认带宽。也许稍微窄一点的带宽会更好一些。（我使用adj=.5了下面的参数。您可以阅读 R 文档density以了解如何更改带宽和内核类型。）

hist(x, prob=T, ylim=c(0, .01), col="skyblue2"); rug(x)
 lines(density(x, adj=.5), col="red")

如果您愿意，您可以获得大约 500 个 KDE 高度的打印输出，然后从左到右扫描数字列表以查找增加和减少以定位模式。（如果必须自动扫描，您可以获取连续高度的差异，并查看它们在何处更改符号以定位模式。）

以下是 的一些输出density，包括前 100 个高度，此处四舍五入到三个位置。

DEN = density(x, adj=.5)
DEN

Call:
        density.default(x = x, adjust = 0.5)

Data: x (950 obs.);     Bandwidth 'bw' = 6.928

       x                y            
 Min.   :-13.61   Min.   :1.154e-06  
 1st Qu.: 54.26   1st Qu.:1.829e-03  
 Median :122.13   Median :3.131e-03  
 Mean   :122.13   Mean   :3.680e-03  
 3rd Qu.:190.00   3rd Qu.:5.655e-03  
 Max.   :257.87   Max.   :1.010e-02  

round(DEN$y[1:100],3)
  [1] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [11] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [21] 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.001 0.001 0.001
 [31] 0.001 0.001 0.001 0.001 0.002 0.002 0.002 0.002 0.002 0.003
 [41] 0.003 0.003 0.004 0.004 0.004 0.005 0.005 0.006 0.006 0.006
 [51] 0.007 0.007 0.007 0.008 0.008 0.009 0.009 0.009 0.009 0.010
 [61] 0.010 0.010 0.010 0.010 0.010 0.010 0.010 0.010 0.010 0.010
 [71] 0.009 0.009 0.009 0.009 0.008 0.008 0.008 0.007 0.007 0.007
 [81] 0.006 0.006 0.006 0.005 0.005 0.005 0.004 0.004 0.004 0.004
 [91] 0.003 0.003 0.003 0.003 0.003 0.002 0.002 0.002 0.002 0.002

不知道你在寻找什么样的答案，但我想我可以试一试。如果您有其中一些，可能最容易手动完成。在最左边的“凹凸”的波谷处切断样品，找到样品模式。

既然你问这个，我假设你有一堆这些并且不能手工完成。这是算法的伪代码：

• 放i = 0
• 从 0 到 的窗口开始i，找出该窗口内数据的模式。
• 如果模式增加，继续增加i并重复上一步。
• 如果模式停止增加，你已经找到了你的“最小”模式。

对于“最大”模式，只需从右侧执行相同操作即可。