2011-06-15 14 views
9

Potrzebuję automatycznie wykrywać spadki na wykresie 2D, podobnie jak obszary oznaczone czerwonymi okręgami na poniższym rysunku. Interesują mnie tylko "główne" spadki, co oznacza, że ​​spadki muszą mieć minimalną długość na osi x. Liczba spadków jest nieznana, tj. Różne wykresy będą zawierały różną liczbę spadków. Jakieś pomysły?Wykrywanie spadków na wykresie 2D

Dips in a 2D plot

Aktualizacja:

Zgodnie z wnioskiem, oto przykładowe dane, wraz z próbą wygładzić za pomocą filtrowania mediana, jak sugeruje winorośli.

Wygląda na to, że potrzebuję teraz solidnego sposobu na przybliżenie pochodnej w każdym punkcie, który zignorowałby małe blipy, które pozostają w danych. Czy istnieje jakieś standardowe podejście?

y <- c(0.9943,0.9917,0.9879,0.9831,0.9553,0.9316,0.9208,0.9119,0.8857,0.7951,0.7605,0.8074,0.7342,0.6374,0.6035,0.5331,0.4781,0.4825,0.4825,0.4879,0.5374,0.4600,0.3668,0.3456,0.4282,0.3578,0.3630,0.3399,0.3578,0.4116,0.3762,0.3668,0.4420,0.4749,0.4556,0.4458,0.5084,0.5043,0.5043,0.5331,0.4781,0.5623,0.6604,0.5900,0.5084,0.5802,0.5802,0.6174,0.6124,0.6374,0.6827,0.6906,0.7034,0.7418,0.7817,0.8311,0.8001,0.7912,0.7912,0.7540,0.7951,0.7817,0.7644,0.7912,0.8311,0.8311,0.7912,0.7688,0.7418,0.7232,0.7147,0.6906,0.6715,0.6681,0.6374,0.6516,0.6650,0.6604,0.6124,0.6334,0.6374,0.5514,0.5514,0.5412,0.5514,0.5374,0.5473,0.4825,0.5084,0.5126,0.5229,0.5126,0.5043,0.4379,0.4781,0.4600,0.4781,0.3806,0.4078,0.3096,0.3263,0.3399,0.3184,0.2820,0.2167,0.2122,0.2080,0.2558,0.2255,0.1921,0.1766,0.1732,0.1205,0.1732,0.0723,0.0701,0.0405,0.0643,0.0771,0.1018,0.0587,0.0884,0.0884,0.1240,0.1088,0.0554,0.0607,0.0441,0.0387,0.0490,0.0478,0.0231,0.0414,0.0297,0.0701,0.0502,0.0567,0.0405,0.0363,0.0464,0.0701,0.0832,0.0991,0.1322,0.1998,0.3146,0.3146,0.3184,0.3578,0.3311,0.3184,0.4203,0.3578,0.3578,0.3578,0.4282,0.5084,0.5802,0.5667,0.5473,0.5514,0.5331,0.4749,0.4037,0.4116,0.4203,0.3184,0.4037,0.4037,0.4282,0.4513,0.4749,0.4116,0.4825,0.4918,0.4879,0.4918,0.4825,0.4245,0.4333,0.4651,0.4879,0.5412,0.5802,0.5126,0.4458,0.5374,0.4600,0.4600,0.4600,0.4600,0.3992,0.4879,0.4282,0.4333,0.3668,0.3005,0.3096,0.3847,0.3939,0.3630,0.3359,0.2292,0.2292,0.2748,0.3399,0.2963,0.2963,0.2385,0.2531,0.1805,0.2531,0.2786,0.3456,0.3399,0.3491,0.4037,0.3885,0.3806,0.2748,0.2700,0.2657,0.2963,0.2865,0.2167,0.2080,0.1844,0.2041,0.1602,0.1416,0.2041,0.1958,0.1018,0.0744,0.0677,0.0909,0.0789,0.0723,0.0660,0.1322,0.1532,0.1060,0.1018,0.1060,0.1150,0.0789,0.1266,0.0965,0.1732,0.1766,0.1766,0.1805,0.2820,0.3096,0.2602,0.2080,0.2333,0.2385,0.2385,0.2432,0.1602,0.2122,0.2385,0.2333,0.2558,0.2432,0.2292,0.2209,0.2483,0.2531,0.2432,0.2432,0.2432,0.2432,0.3053,0.3630,0.3578,0.3630,0.3668,0.3263,0.3992,0.4037,0.4556,0.4703,0.5173,0.6219,0.6412,0.7275,0.6984,0.6756,0.7079,0.7192,0.7342,0.7458,0.7501,0.7540,0.7605,0.7605,0.7342,0.7912,0.7951,0.8036,0.8074,0.8074,0.8118,0.7951,0.8118,0.8242,0.8488,0.8650,0.8488,0.8311,0.8424,0.7912,0.7951,0.8001,0.8001,0.7458,0.7192,0.6984,0.6412,0.6516,0.5900,0.5802,0.5802,0.5762,0.5623,0.5374,0.4556,0.4556,0.4333,0.3762,0.3456,0.4037,0.3311,0.3263,0.3311,0.3717,0.3762,0.3717,0.3668,0.3491,0.4203,0.4037,0.4149,0.4037,0.3992,0.4078,0.4651,0.4967,0.5229,0.5802,0.5802,0.5846,0.6293,0.6412,0.6374,0.6604,0.7317,0.7034,0.7573,0.7573,0.7573,0.7772,0.7605,0.8036,0.7951,0.7817,0.7869,0.7724,0.7869,0.7869,0.7951,0.7644,0.7912,0.7275,0.7342,0.7275,0.6984,0.7342,0.7605,0.7418,0.7418,0.7275,0.7573,0.7724,0.8118,0.8521,0.8823,0.8984,0.9119,0.9316,0.9512) 

yy <- runmed(y, 41) 
plot(y, type="l", ylim=c(0,1), ylab="", xlab="", lwd=0.5) 
points(yy, col="blue", type="l", lwd=2) 

Median filtering

+0

Chyba można wygładzić dane trochę i użyj tego: http://stackoverflow.com/questions/6324354/ add-a-curve-that-fits-the-peaks-from-a-plot-in-r/ –

+2

dane przykładowe byłyby miłe ... –

+1

@Joris Dodałem dane użyte do wygenerowania wykresu. Dziękuję za wskazanie tego. – Leo

Odpowiedz

6

EDYCJA: Funkcja usuwa regiony, aby zawierały tylko najniższą część, jeśli jest to potrzebne.

W rzeczywistości użycie średniej jest łatwiejsze niż użycie mediany. Dzięki temu można znaleźć regiony, w których rzeczywiste wartości są stale poniżej średniej. Mediana nie jest wystarczająco gładka, aby można ją było łatwo zastosować.

Przykładem funkcji, aby zrobić to byłoby:

FindLowRegion <- function(x,n=length(x)/4,tol=length(x)/20,p=0.5){ 
    nx <- length(x) 
    n <- 2*(n %/% 2) + 1 
    # smooth out based on means 
    sx <- rowMeans(embed(c(rep(NA,n/2),x,rep(NA,n/2)),n),na.rm=T) 
    # find which series are far from the mean 
    rlesx <- rle((sx-x)>0) 
    # construct start and end of regions 
    int <- embed(cumsum(c(1,rlesx$lengths)),2) 
    # which regions fulfill requirements 
    id <- rlesx$value & rlesx$length > tol 
    # Cut regions to be in general smaller than median 
    regions <- 
    apply(int[id,],1,function(i){ 
     i <- min(i):max(i) 
     tmp <- x[i] 
     id <- which(tmp < quantile(tmp,p)) 
     id <- min(id):max(id) 
     i[id]    
    }) 
    # return 
    unlist(regions) 
} 

gdzie

  • n określa ile wartości są wykorzystywane do obliczania biegu myśli,
  • tol określa, ile kolejnych wartości powinny być niższym niż średni bieg, aby mówić o niskim regionie, i
  • p określa zastosowany odcięcie (jako kwantyl) do rozebrania regionów do ich najniższej części. Kiedy p = 1, pokazany jest kompletny dolny region.

Funkcja została zmodyfikowana tak, aby działała zgodnie z prezentowanymi danymi, ale liczby mogą wymagać niewielkiej korekty w celu pracy z innymi danymi.

Ta funkcja zwraca zestaw indeksów, który pozwala znaleźć najniższe regiony.Ilustrowane swojej y wektora:

Lows <- FindLowRegion(y) 

newx <- seq_along(y) 
newy <- ifelse(newx %in% Lows,y,NA) 
plot(y, col="blue", type="l", lwd=2) 
lines(newx,newy,col="red",lwd="3") 

Daje:

enter image description here

+0

Czy ważność statystyczna? – hadley

+1

@Hadley: o tyle samo, co less() i przyjaciele. –

+1

W twoim stwierdzeniu 'newy <- ifelse (x% w% Lows, y, NA)', skąd pochodzi 'x'? Czy nie powinno to być "newx"? –

3

trzeba wygładzić wykres w jakiś sposób. Median filtration jest bardzo przydatny do tego celu (patrz http://en.wikipedia.org/wiki/Median_filter). Po wygładzeniu będziesz musiał po prostu szukać minimów, tak jak zwykle (tj. Szukać punktów, w których pierwsza pochodna zmienia się z ujemnej na pozytywną).

+0

Dzięki za sugestię. Zaktualizowałem pytanie przy pomocy filtrowania mediany. Nadal jest trochę hałasu, więc pozostaje jeden problem: dokładne przybliżenie 1. pochodnej. – Leo

+0

@Leo: Znam praktycznie nic o R ... Ale mówiąc algorytmicznie, pierwszą rzeczą, którą chciałbym spróbować jest okno przesuwne: gdy wszystkie punkty w oknie, z wyjątkiem skrajnych lewych i prawych, są poniżej najniższego z lewej i rightmost, a następnie dip jest znaleziony, a okno jest przesuwane o jego szerokość na raz, w przeciwnym razie okno jest przesuwane o jeden krok ... – vines

+1

Średnie wygładzenie jest o wiele bardziej przydatne, jak pokazano w mojej odpowiedzi. –

0

Moja pierwsza myśl była czymś znacznie okrutniejszym niż filtrowanie. Dlaczego nie szukać wielkich kropli, a następnie wystarczająco długich, stabilnych okresów?

span.b <- 20 
threshold.b <- 0.2 
dy.b <- c(rep(NA, span.b), diff(y, lag = span.b)) 
span.f <- 10 
threshold.f <- 0.05 
dy.f <- c(diff(y, lag = span.f), rep(NA, span.f)) 
down <- which(dy.b < -1 * threshold.b & abs(dy.f) < threshold.f) 
abline(v = down) 

Wykres pokazuje, że nie jest doskonały, ale nie odrzucić odstających (myślę, że to zależy od Wciel się w danych).

1

Prostsza odpowiedź (która również nie wymaga wygładzania) mogą być dostarczone poprzez dostosowanie funkcji maxdrawdown() z tseries. A drawdown jest powszechnie definiowany jako odwrót od najnowszego maksimum; tutaj chcemy czegoś przeciwnego. Taką funkcję można następnie wykorzystać w przesuwającym oknie nad danymi lub nad danymi posegmentowanymi.

maxdrawdown <- function(x) { 
    if(NCOL(x) > 1) 
     stop("x is not a vector or univariate time series") 
    if(any(is.na(x))) 
     stop("NAs in x") 
    cmaxx <- cummax(x)-x 
    mdd <- max(cmaxx) 
    to <- which(mdd == cmaxx) 
    from <- double(NROW(to)) 
    for (i in 1:NROW(to)) 
     from[i] <- max(which(cmaxx[1:to[i]] == 0)) 
    return(list(maxdrawdown = mdd, from = from, to = to)) 
} 

Więc zamiast korzystania cummax(), trzeba by przejść do cummin() itp

+0

Lubię proste odpowiedzi. – Jubbles