寻找局部最大值和最小值

Question

我正在寻找一种计算有效的方法来为R中的大数字列表找到局部最大值/最小值.希望没有for循环...

例如,如果我有一个数据文件1 2 3 2 1 1 2 1,我希望函数返回3和7,它们是局部最大值的位置.

Answer 1

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/)适用于数据较嘈杂的情况.

Answer 2

@ 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

Answer 3

使用动物园库函数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软件包(用于进行质谱分析的"峰值概率对比"),因为在阅读@BenBolker上面的评论之前我没有意识到它的可用性,而且我认为添加这几个字会增加有人的机会.大规模的兴趣将在搜索中看到这一点.)

Answer 4

提供了一些很好的解决方案,但这取决于您的需求.

只是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
         }
     }
 }

Answer 5

我今天抓了这个.我知道你希望没有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)

Answer 6

聚会迟到了，但其他人可能会对此感兴趣。你可以使用目前的（内部）函数find_peaks从ggpmisc包。您可以使用threshold,span和strictarguments对其进行参数化。由于ggpmisc包旨在与ggplot2您一起使用，您可以使用和函数直接绘制最小值和最大值：stat_peaksstat_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")

Answer 7

@42- 的回答很棒，但我有一个用例，我不想使用zoo. dplyr使用lagand很容易实现这一点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解决方案一样，您可以分别通过lag/lead参数n和控制窗口大小和边缘情况default。