我正在寻找一种在R中对大量数字进行局部最大值/最小值计算的高效算法。
希望不使用
for循环...
例如，如果我有一个数据文件如1 2 3 2 1 1 2 1，我希望该函数返回3和7，它们是局部最大值的位置。17个回答
70
diff(diff(x))（或者diff(x,differences=2)：感谢@ZheyuanLi）本质上计算了二阶导数的离散模拟，因此在局部极大值处应该为负。下面的+1用于处理diff结果长度小于输入向量的情况。
编辑：添加了@Tommy关于delta-x不等于1的修正...
tt <- c(1,2,3,2,1, 1, 2, 1)
which(diff(sign(diff(tt)))==-2)+1
我的建议（http://statweb.stanford.edu/~tibs/PPC/Rdist/）适用于数据较嘈杂的情况。
- Ben Bolker
7
44
@Ben的解决方案非常棒。不过它无法处理以下情况：
# all these return numeric(0):
x <- c(1,2,9,9,2,1,1,5,5,1) # duplicated points at maxima
which(diff(sign(diff(x)))==-2)+1
x <- c(2,2,9,9,2,1,1,5,5,1) # duplicated points at start
which(diff(sign(diff(x)))==-2)+1
x <- c(3,2,9,9,2,1,1,5,5,1) # start is maxima
which(diff(sign(diff(x)))==-2)+1
这里是一个更健壮（但更慢，更丑）的版本：
localMaxima <- function(x) {
# Use -Inf instead if x is numeric (non-integer)
y <- diff(c(-.Machine$integer.max, x)) > 0L
rle(y)$lengths
y <- cumsum(rle(y)$lengths)
y <- y[seq.int(1L, length(y), 2L)]
if (x[[1]] == x[[2]]) {
y <- y[-1]
}
y
}
x <- c(1,2,9,9,2,1,1,5,5,1)
localMaxima(x) # 3, 8
x <- c(2,2,9,9,2,1,1,5,5,1)
localMaxima(x) # 3, 8
x <- c(3,2,9,9,2,1,1,5,5,1)
localMaxima(x) # 1, 3, 8
- Tommy
5
谢谢，我尝试了这段代码，它可以工作！你如何在不改变输入的情况下修改它以适应局部最小值？ - Vahid Mirjalili
嗨Tommy，我想在一个包中使用你的localMinima函数，你能联系我一下吗，这样我就可以适当地致谢了吗？ - Yann Abraham
@VahidMir 基本上，这个函数是一种（聪明！）的方法，用于获取向量的一阶导数从正变为负的位置。因此，局部最小值将在从负变为正的位置给出：只需将第一行替换为
y <- diff(c(.Machine$integer.max, x)) < 0L（这保留了检测初始最小值的可能性）。 - ztl4
localMaxima() 函数在拐点处产生误报，例如 localMaxima(c(1, 2, 2, 3, 2, 1)) 返回的结果是 2 4 而不是仅返回 4。需要进行改进。 - jacanterbury为什么要两次调用rle(y)$lengths？我的意思是，我理解
y <- cumsum(rle(y)$lengths)，但不理解前面的独立的rle(y)$lengths。 - Andy27
使用zoo库的rollapply函数：
x <- c(1, 2, 3, 2, 1, 1, 2, 1)
library(zoo)
xz <- as.zoo(x)
rollapply(xz, 3, function(x) which.min(x)==2)
# 2 3 4 5 6 7
#FALSE FALSE FALSE TRUE FALSE FALSE
rollapply(xz, 3, function(x) which.max(x)==2)
# 2 3 4 5 6 7
#FALSE TRUE FALSE FALSE FALSE TRUE
使用“coredata”提取那些'which.max'是表示局部最大值的“中心值”的值的索引。您可以显然使用which.min而不是which.max来查找局部最小值。
rxz <- rollapply(xz, 3, function(x) which.max(x)==2)
index(rxz)[coredata(rxz)]
#[1] 3 7
我假设你不想要起始值或结束值，但如果你需要的话，你可以在处理之前填充向量的末尾，就像端粒在染色体上的作用一样。
（我注意到了ppc包（“Peak Probability Contrasts”）用于进行质谱分析，只是因为在阅读@BenBolker上面的评论之前，我不知道它的可用性，并且我认为增加这几个词将增加质谱学感兴趣的人通过搜索找到这个问题。）
- IRTFM
7
3相对于其他方法，这种方法有非常明显的优势。将间隔增加到大于3的某个值后，即使存在其他更接近的点，我们也可以忽略某个点略微高于其两个最近邻居的情况。这在测量数据具有小的随机变化时非常有用。 - jpmc26
2这是一个很棒的解决方案，但需要注意的是：最好明确定义
align 参数。zoo:::rollapply.zoo 默认使用 align = "center"，但 xts:::rollapply.xts 使用 align = "right"。 - mikeck5@dleal，你需要在数组
xz上滚动一个宽度为3的窗口。这个窗口的内容是函数参数x，该函数返回最大值的索引。如果该索引指向窗口中心，则表示你已经找到了局部最大值！在这种特殊情况下，窗口宽度为3，因此中间元素的索引为2。基本上，你需要寻找一个条件which.max(x) == m，其中窗口宽度等于2*m-1。 - R Kiselev1有趣的是，@42-的建议在您的重复值比宽度（例如3）更多的情况下会失败。换句话说，一个鞍点会被误认为是一个极值。一个简单的例子是：
x <- c(3, 2, 2, 2, 2, 1, 3)，然后 rx <- rollapply(as.zoo(x), 3, function(x) {which.min(x)==2})，而 index(rx)[coredata(rx)] 错误地给出了 [1] 2 6（正确应该是 [1] 6）。 - user3375672使用间隔为5（
rxz <- rollapply(xz, 5, function(x) which.max(x)==2)）添加数据到@jpmc26，例如x <- c(1, 2, 3, 2, 3, 4, 5, 6, 6, 5 ,4, 3, 2)，早期的噪声干扰3,2,3被忽略，但高原6,6被正确检测。 - jacanterbury显示剩余2条评论
19
今天我试着做了一下这个。我知道你说尽量不要用for循环，但我坚持使用了apply函数。它相对紧凑快速，并允许阈值规定，所以你可以大于1。
该函数：
inflect <- function(x, threshold = 1){
up <- sapply(1:threshold, function(n) c(x[-(seq(n))], rep(NA, n)))
down <- sapply(-1:-threshold, function(n) c(rep(NA,abs(n)), x[-seq(length(x), length(x) - abs(n) + 1)]))
a <- cbind(x,up,down)
list(minima = which(apply(a, 1, min) == a[,1]), maxima = which(apply(a, 1, max) == a[,1]))
}
要可视化/玩弄阈值，您可以运行以下代码：
# Pick a desired threshold # to plot up to
n <- 2
# Generate Data
randomwalk <- 100 + cumsum(rnorm(50, 0.2, 1)) # climbs upwards most of the time
bottoms <- lapply(1:n, function(x) inflect(randomwalk, threshold = x)$minima)
tops <- lapply(1:n, function(x) inflect(randomwalk, threshold = x)$maxima)
# Color functions
cf.1 <- grDevices::colorRampPalette(c("pink","red"))
cf.2 <- grDevices::colorRampPalette(c("cyan","blue"))
plot(randomwalk, type = 'l', main = "Minima & Maxima\nVariable Thresholds")
for(i in 1:n){
points(bottoms[[i]], randomwalk[bottoms[[i]]], pch = 16, col = cf.1(n)[i], cex = i/1.5)
}
for(i in 1:n){
points(tops[[i]], randomwalk[tops[[i]]], pch = 16, col = cf.2(n)[i], cex = i/1.5)
}
legend("topleft", legend = c("Minima",1:n,"Maxima",1:n),
pch = rep(c(NA, rep(16,n)), 2), col = c(1, cf.1(n),1, cf.2(n)),
pt.cex = c(rep(c(1, c(1:n) / 1.5), 2)), cex = .75, ncol = 2)
- Evan Friedland
4
我喜欢函数<inflect>。感谢Evan (y)。 - Hoang Le
很高兴听到这个消息！ - Evan Friedland
我的数据存在最大值和最小值的重复值，例如c(1,2,3,3,2,1,1,2,3)，这会为每个最大/最小值生成多个命中。尝试使用“阈值”进行实验似乎只会改变图表上的点大小，但无法解决这个问题。有什么建议吗？ - jacanterbury
你好，当数值相等时，你希望发生什么？该函数将它们视为相同的点--应该被忽略。在一个平坦的山顶上，无论你走到哪里都是山峰，即使中间有裂缝。你正在处理什么样的数据，使得它们的值相等？另外，如果相邻高度不相同，点的变化也会不同。请参见向量
c(0,0,0,1,0.7,3,2,3,3,2,1,1,2,3,0.7, 0.5,0,0,0)，阈值为3。 - Evan Friedland13
有一些不错的解决方案，但这取决于您的需求。
仅使用diff(tt)即可返回差异。
您希望检测何时从递增值转变为递减值。其中一种方法由@Ben提供：
diff(sign(diff(tt)))==-2
这里的问题是它只能检测到从严格增加到严格减少的变化。
稍作更改即可允许峰值处出现重复值（返回峰值最后一次出现的TRUE）：
diff(diff(x)>=0)<0
如果您想检测开头或结尾的极大值，那么您只需要正确地填充前面和后面即可。
以下是包括寻找谷底在内的所有内容的函数：
which.peaks <- function(x,partial=TRUE,decreasing=FALSE){
if (decreasing){
if (partial){
which(diff(c(FALSE,diff(x)>0,TRUE))>0)
}else {
which(diff(diff(x)>0)>0)+1
}
}else {
if (partial){
which(diff(c(TRUE,diff(x)>=0,FALSE))<0)
}else {
which(diff(diff(x)>=0)<0)+1
}
}
}
- jamie.f.olson
1
1亲爱的未来访客，我尝试了这里提供的几个解决方案，其中这个对我来说效果最好。 - jacanterbury
13
虽然有些晚，但这对其他人可能很有兴趣。现在你可以使用ggpmisc软件包中的（内部）函数find_peaks。您可以使用threshold、span和strict参数对其进行参数化。由于ggpmisc软件包旨在与ggplot2一起使用，因此您可以直接使用stat_peaks和stat_valleys函数绘制极小值和极大值：
set.seed(1)
x <- 1:10
y <- runif(10)
# Maxima
x[ggpmisc:::find_peaks(y)]
[1] 4 7
y[ggpmisc:::find_peaks(y)]
[1] 0.9082078 0.9446753
# Minima
x[ggpmisc:::find_peaks(-y)]
[1] 5
y[ggpmisc:::find_peaks(-y)]
[1] 0.2016819
# Plot
ggplot(data = data.frame(x, y), aes(x = x, y = y)) + geom_line() + stat_peaks(col = "red") + stat_valleys(col = "green")
- symbolrush
2
谢谢你的回复。但是我如何在ggplot2环境之外使用stat_peaks呢？ - undefined
在我看来，如果在
ggplot2环境之外使用stat_peaks并没有太多意义。因为stat_peaks旨在在ggplot调用内使用。另一方面，我回答的第一部分在任何地方都适用。使用ggpmisc:::find_peaks计算极小值和极大值是通用的。 - undefined8
@42-的回答很好，但是我有一个使用案例，不想使用zoo。使用dplyr和lag、lead很容易实现：
library(dplyr)
test = data_frame(x = sample(1:10, 20, replace = TRUE))
mutate(test, local.minima = if_else(lag(x) > x & lead(x) > x, TRUE, FALSE)
与 rollapply 解决方案类似，您可以通过 n 和 default 参数来控制窗口大小和边缘情况的 lag/lead 方法。
- mikeck
2
当窗口大于1时，此解决方案与rollapply解决方案之间存在差异。假设我们想要在任一方向上查看3个位置。使用rollapply解决方案，我们可以查看7个值，并且它会告诉我们中间的值是否为最小值。在此解决方案中，使用if_else（lag（x，3）> x＆lead（x，3）> x）只会查看第三个位置，而不是1和2。我喜欢使用dplyr解决方案的想法，但编写6个条件似乎有点繁琐。 - canderson156
请看我的答案，它是在此基础上构建的解决方案。 - canderson156
4
在我所处理的情况中，重复数据经常出现。因此，我已经实现了一个函数，允许查找第一个或最后一个极值（最小值或最大值）：
locate_xtrem <- function (x, last = FALSE)
{
# use rle to deal with duplicates
x_rle <- rle(x)
# force the first value to be identified as an extrema
first_value <- x_rle$values[1] - x_rle$values[2]
# differentiate the series, keep only the sign, and use 'rle' function to
# locate increase or decrease concerning multiple successive values.
# The result values is a series of (only) -1 and 1.
#
# ! NOTE: with this method, last value will be considered as an extrema
diff_sign_rle <- c(first_value, diff(x_rle$values)) %>% sign() %>% rle()
# this vector will be used to get the initial positions
diff_idx <- cumsum(diff_sign_rle$lengths)
# find min and max
diff_min <- diff_idx[diff_sign_rle$values < 0]
diff_max <- diff_idx[diff_sign_rle$values > 0]
# get the min and max indexes in the original series
x_idx <- cumsum(x_rle$lengths)
if (last) {
min <- x_idx[diff_min]
max <- x_idx[diff_max]
} else {
min <- x_idx[diff_min] - x_rle$lengths[diff_min] + 1
max <- x_idx[diff_max] - x_rle$lengths[diff_max] + 1
}
# just get number of occurences
min_nb <- x_rle$lengths[diff_min]
max_nb <- x_rle$lengths[diff_max]
# format the result as a tibble
bind_rows(
tibble(Idx = min, Values = x[min], NB = min_nb, Status = "min"),
tibble(Idx = max, Values = x[max], NB = max_nb, Status = "max")) %>%
arrange(.data$Idx) %>%
mutate(Last = last) %>%
mutate_at(vars(.data$Idx, .data$NB), as.integer)
}
对于原始问题的答案是:
> x <- c(1, 2, 3, 2, 1, 1, 2, 1)
> locate_xtrem(x)
# A tibble: 5 x 5
Idx Values NB Status Last
<int> <dbl> <int> <chr> <lgl>
1 1 1 1 min FALSE
2 3 3 1 max FALSE
3 5 1 2 min FALSE
4 7 2 1 max FALSE
5 8 1 1 min FALSE
结果表明第二个最小值等于1并且该值从索引5开始重复两次。因此，通过指示函数查找局部极值的最后出现可以获得不同的结果：
> locate_xtrem(x, last = TRUE)
# A tibble: 5 x 5
Idx Values NB Status Last
<int> <dbl> <int> <chr> <lgl>
1 1 1 1 min TRUE
2 3 3 1 max TRUE
3 6 1 2 min TRUE
4 7 2 1 max TRUE
5 8 1 1 min TRUE
根据目标，可以在局部极值的第一个和最后一个值之间进行切换。使用last=TRUE得到的第二个结果也可以通过"Idx"和"NB"两列之间的运算得到...
最后，为了处理数据中的噪声，可以实现一个函数来删除低于给定阈值的波动。由于这超出了初始问题的范围，因此不公开代码。我已将其封装在一个包中（主要是为了自动化测试过程），以下是一个结果示例：
x_series %>% xtrem::locate_xtrem()
x_series %>% xtrem::locate_xtrem() %>% remove_noise()
- Jeffery Petit
2
不错！我有一个累积值图表，因此有平坦的范围。你的解决方案是唯一有效的。如果能包含绘图代码就更好了。 - James Hirschorn
这是一个很棒的函数！能够在重复序列中控制首尾非常方便。 - vb66
2
之前的解决方案中，我在获取位置信息时遇到了一些问题，并想出了一种直接获取最小值和最大值的方法。以下代码将实现此功能，并绘制图表，在图表中用绿色标记最小值，用红色标记最大值。与which.max()函数不同的是，这将从数据框中提取所有最小值/最大值的索引。第一个diff()函数中添加零值是为了解决使用该函数时会导致结果长度减少的问题。将其插入到最内层的diff()函数调用中可以避免在逻辑表达式外部添加偏移量。虽然这并不重要，但我觉得这是一种更清晰的方法。
# create example data called stockData
stockData = data.frame(x = 1:30, y=rnorm(30,7))
# get the location of the minima/maxima. note the added zero offsets
# the location to get the correct indices
min_indexes = which(diff( sign(diff( c(0,stockData$y)))) == 2)
max_indexes = which(diff( sign(diff( c(0,stockData$y)))) == -2)
# get the actual values where the minima/maxima are located
min_locs = stockData[min_indexes,]
max_locs = stockData[max_indexes,]
# plot the data and mark minima with red and maxima with green
plot(stockData$y, type="l")
points( min_locs, col="red", pch=19, cex=1 )
points( max_locs, col="green", pch=19, cex=1 )
- Ehren
4
几乎非常好 - 似乎无法在末尾使用最大值> histData $ counts
[1] 18000 0 0 0 0 0 0 0 0 0 0 0
[217] 0 0 0 0 0 0 0 0 5992 - idontgetoutmuch
max_indexes = sign(diff( c(0,histData$counts,0)))) 这个代码虽然能正常工作，但我不知道它是否会影响其他部分。 - idontgetoutmuch@idontgetoutmuch... 这个方法基本上使用数据的一阶导数计算，不会在被评估的系列的端点找到相对最大值或最小值。如果倒数第二个值是相对最大/最小值，那么它将适用，因为可以在那里近似导数。如果你正在寻找系列中的最大值，max()函数应该可以很好地工作。结合上面的代码，应该能够获得有关最大值/最小值的信息。 - Ehren
我在上面的评论中第一句话应该更清楚一些... 该方法基本上使用数据的一阶导数逼近，并且不会在被评估的系列的端点处找到相对极大值或极小值，因为无法知道端点是否是相对极大/小值。 - Ehren
2
这里是minima主题的解决方案：
@Ben的解决方案
请参考Tommy的帖子中的案例！
@Tommy的解决方案：
请注意：无论是
@Ben的解决方案
x <- c(1,2,3,2,1,2,1)
which(diff(sign(diff(x)))==+2)+1 # 5
请参考Tommy的帖子中的案例！
@Tommy的解决方案：
localMinima <- function(x) {
# Use -Inf instead if x is numeric (non-integer)
y <- diff(c(.Machine$integer.max, x)) > 0L
rle(y)$lengths
y <- cumsum(rle(y)$lengths)
y <- y[seq.int(1L, length(y), 2L)]
if (x[[1]] == x[[2]]) {
y <- y[-1]
}
y
}
x <- c(1,2,9,9,2,1,1,5,5,1)
localMinima(x) # 1, 7, 10
x <- c(2,2,9,9,2,1,1,5,5,1)
localMinima(x) # 7, 10
x <- c(3,2,9,9,2,1,1,5,5,1)
localMinima(x) # 2, 7, 10
请注意：无论是
localMaxima还是localMinima，都无法处理起始位置的重复最大值/最小值！- Sebastian
6
不确定你的答案真正能为问题带来什么，因为其他答案已经包含了相同的算法。 - m4rtin
2是的，但它提供了最小值的解决方案，就像最初所问的那样。此外，还没有提到起始点处重复的最大值/最小值情况。 - Sebastian
好吧...即使这基本上是相同的答案，我也无法反驳。所以我不会点踩，但也不会点赞。你应该尝试回答那些还没有被回答的问题（即使这个问题还没有正式回答，但前两个答案已经有18和20票了，这就是一样的）。 - m4rtin
顺便说一句，我可能需要帮助找到一种适应更大间隔的最大/最小值函数的方法，因此“delta x > 1”。有人有想法吗？ - Sebastian
@Sebastian，如果你还没有看到的话，请看一下我关于更大间隔的答案。 - Evan Friedland
@Evan 确实是一个不错的解决方案，感谢您还提供了漂亮的图示！ - Sebastian
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接
which(diff(sign(diff(x)))==-2)+1。 - Tommyabs()适当地，如果你想要同时寻找两者），你应该将==-2更改为==2。 - Ben Bolker