Я ищу эффективный вычислительно эффективный способ найти локальные максимумы/минимумы для большого списка чисел в R.
Надеемся, что без петель for...
Например, если у меня есть файл данных типа 1 2 3 2 1 1 2 1, я хочу, чтобы функция возвращала 3 и 7, которые являются положениями локальных максимумов.
diff(diff(x)) (или diff(x,differences=2): благодаря @ZheyuanLi) по существу вычисляет дискретный аналог второй производной, поэтому должен быть отрицательным при локальных максимумах. Ниже +1 заботится о том, что результат diff короче входного вектора.
edit: добавлена коррекция @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 решение довольно сладкое. Однако он не обрабатывает следующие случаи:
# 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
Использовать функцию библиотеки zoo:
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 выше, и я думаю, что добавление этих нескольких слов увеличит шансы, что кто-то с массовым интересом увидит это при поиске.)
Сегодня я сделал удар. Я знаю, что вы сказали, что без петель, но я застрял в использовании функции 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)
Есть несколько хороших решений, но это зависит от того, что вам нужно.
Просто 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
}
}
}
Ответа на этот вопрос @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, вы можете управлять размерами окна и краями через аргументы lag/lead n и default соответственно.
Здесь решение для минимумов:
@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 не могут обрабатывать дублированные максимумы/минимумы при запуске!
У меня были проблемы с работой в предыдущих решениях, и я попытался получить максимальные минимумы и максимумы. Приведенный ниже код сделает это и закроет его, отметив минимальные зеленые и максимумы красным. В отличие от функции 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 )
В пакете pracma используйте
tt <- c(1,2,3,2,1, 1, 2, 1)
tt_peaks <- findpeaks(tt, zero = "0", peakpat = NULL,
minpeakheight = -Inf, minpeakdistance = 1, threshold = 0, npeaks = 0, sortstr = FALSE)
[,1] [,2] [,3] [,4]
[1,] 3 3 1 5
[2,] 2 7 6 8
Это возвращает матрицу с 4 столбцами. Первый столбец показывает абсолютные значения локальных пиков. Второй столбец - это индексы. Третий и четвертый столбцы - это начало и конец пиков (с потенциальным перекрытием).
Подробнее см. Https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/findpeaks.
Одно предостережение: я использовал его в серии нецелых чисел, и пик был одним индексом слишком поздним (для всех пиков), и я не знаю почему. Поэтому мне пришлось вручную удалить "1" из моего индексного вектора (неважно).
Я разместил это в другом месте, но я думаю, что это интересный способ сделать это. Я не уверен, какова его вычислительная эффективность, но это очень сжатый способ решения проблемы.
vals=rbinom(1000,20,0.5)
text=paste0(substr(format(diff(vals),scientific=TRUE),1,1),collapse="")
sort(na.omit(c(gregexpr('[ ]-',text)[[1]]+1,ifelse(grepl('^-',text),1,NA),
ifelse(grepl('[^-]$',text),length(vals),NA))))