국소 최댓값과 최솟값 찾기
R에서 큰 숫자 목록에 대한 로컬 최대 값 / 최소값을 찾는 계산적으로 효율적인 방법을 찾고 for있습니다. 루프가 없으면 좋겠습니다.
예를 들어,와 같은 데이터 파일이있는 경우 1 2 3 2 1 1 2 1함수가 로컬 최대 값의 위치 인 3과 7을 반환하기를 원합니다.
diff(diff(x))(또는 diff(x,differences=2): @ZheyuanLi 덕분에) 본질적으로 2 차 미분의 이산 아날로그를 계산하므로 로컬 최대 값에서 음수가되어야합니다. 는 +1아래의 결과가 사실을 담당 diff입력 벡터보다 짧다.
편집 : delta-x가 1이 아닌 경우 @Tommy의 수정을 추가했습니다.
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
동물원 라이브러리 함수 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
그런 다음 'which.max'가 로컬 최대 값을 나타내는 "중심 값"인 값에 대해 'coredata'를 사용하여 인덱스를 가져옵니다. 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의 의견을 읽을 때까지 가용성을 알지 못했기 때문이며,이 몇 단어를 추가하면 a 대량 스펙 관심은 검색에서 볼 수 있습니다.)
몇 가지 좋은 솔루션이 제공되지만 필요한 사항에 따라 다릅니다.
그냥 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
}
}
}
나는 오늘 이것을 찔렀다. 나는 당신이 for 루프없이 희망적으로 말했지만 적용 기능을 사용하는 것을 고수했습니다. 다소 간결하고 빠르며 임계 값 사양을 허용하므로 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)
@ 42-의 답변은 훌륭하지만 사용하고 싶지 않은 사용 사례가있었습니다 zoo. 및 dplyr사용하여 쉽게 구현할 수 있습니다 .laglead
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()함수를 사용할 때마다 발생하는 결과의 누락 된 감소 길이를 설명하기 위해 첫 번째 함수 에 0 값이 추가됩니다 . 이것을 가장 안쪽의 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 )
나는 이것을 다른 곳에 게시했지만 이것은 그것에 대해 흥미로운 방법이라고 생각합니다. 계산 효율성이 무엇인지 잘 모르겠지만 문제를 해결하는 매우 간결한 방법입니다.
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))))
에서 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"을 수동으로 제거해야했습니다.
파티에 늦었지만 이것은 다른 사람들에게 흥미로울 수 있습니다. 요즘에는 패키지 의 (내부) 기능 find_peaks을 사용할 수 있습니다 ggpmisc. 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")
쉽지 않은 시퀀스에 대한 로컬 최대 값과 최소값 찾기 예를 들어 최대 값1 0 1 1 2 0 1 1 0 1 1 1 0 1 은 (1), 5, 7.5, 11 및 (14), 최소값은 2, 6, 9, 13에 위치를 제공합니다.
#Position 1 1 1 1 1
# 1 2 3 4 5 6 7 8 9 0 1 2 3 4
x <- c(1,0,1,1,2,0,1,1,0,1,1,1,0,1) #Frequency
# p v p v p v p v p p..Peak, v..Valey
peakPosition <- function(x, inclBorders=TRUE) {
if(inclBorders) {y <- c(min(x), x, min(x))
} else {y <- c(x[1], x)}
y <- data.frame(x=sign(diff(y)), i=1:(length(y)-1))
y <- y[y$x!=0,]
idx <- diff(y$x)<0
(y$i[c(idx,F)] + y$i[c(F,idx)] - 1)/2
}
#Find Peaks
peakPosition(x)
#1.0 5.0 7.5 11.0 14.0
#Find Valeys
peakPosition(-x)
#2 6 9 13
peakPosition(c(1,2,3,2,1,1,2,1)) #3 7
참고 URL : https://stackoverflow.com/questions/6836409/finding-local-maxima-and-minima
'program story' 카테고리의 다른 글
| Python의 스레드 로컬 저장소 (0) | 2020.11.30 |
|---|---|
| 자바 스크립트에서 날짜가 같은지 확인 (0) | 2020.11.30 |
| 내 로컬 Git 저장소를 원격 Git 저장소로 이동하려면 어떻게해야합니까? (0) | 2020.11.30 |
| ALTER TABLE ADD COLUMN은 시간이 오래 걸립니다. (0) | 2020.11.29 |
| FrameLayout은 무엇을합니까? (0) | 2020.11.29 |