Znajdowanie ekstremów lokalnych funkcji gęstości za pomocą splajnów

15

Próbuję znaleźć lokalne maksima dla funkcji gęstości prawdopodobieństwa (znalezionej metodą R density). Nie mogę wykonać prostej metody „rozglądania się po sąsiadach” (gdzie ktoś rozgląda się po punkcie, aby sprawdzić, czy jest to maksimum lokalne w stosunku do swoich sąsiadów), ponieważ istnieje duża ilość danych. Co więcej, wydaje się bardziej wydajne i ogólne użycie czegoś takiego jak interpolacja splajnu, a następnie znalezienie pierwiastków pierwszej pochodnej, w przeciwieństwie do budowania „rozejrzenia się po sąsiadach” z tolerancją na awarie i innymi parametrami.

Więc moje pytania:

  1. Biorąc pod uwagę funkcję splinefun, jakie metody znajdą lokalne maksima?
  2. Czy istnieje prosty / standardowy sposób na znalezienie pochodnych funkcji zwracanych za pomocą splinefun?
  3. Czy istnieje lepszy / standardowy sposób na znalezienie lokalnych maksimów funkcji gęstości prawdopodobieństwa?

Dla porównania poniżej znajduje się wykres mojej funkcji gęstości. Inne funkcje gęstości, z którymi pracuję, mają podobną formę. Powinienem powiedzieć, że jestem nowy w R, ale nie nowy w programowaniu, więc może istnieć standardowa biblioteka lub pakiet do osiągnięcia tego, czego potrzebuję. funkcja gęstości

Dzięki za pomoc!!

r pdf splines maximum aaronlevin
źródło
Nie jestem pewien, dlaczego duża ilość danych stanowi problem dla metody „rozejrzyj się po sąsiadach”. density()nie szacuje gęstości dla każdego układu odniesienia, szacuje gęstość przy wartościach n , gdzie n jest parametrem określonym przez użytkownika o wartości domyślnej n = 512.
onestop
Moje n to 2 ^ 15 i wydaje się, że dane mają dużą wariancję na poziomie punkt po punkcie. Próbowałem napisać wyszukiwarkę maks./min przy użyciu czegoś podobnego do metody sąsiedztwa (via msExtrema {msProcess}) i byłem w stanie zidentyfikować tylko niektóre maksima, nigdy nie wszystkie, grając z ustawieniami tolerancji.
aaronlevin,
2
Patrząc na kod msExtrema, jest to proste opakowanie peaksz splus2Rpakietu, którego lepiej byłoby użyć bezpośrednio, jeśli chcesz tylko lokalne maksima, a nie lokalne minima. Nie rozumiem, dlaczego użycie wartości domyślnej span=3nie znalazłoby wszystkich lokalnych maksimów. A 2 ^ 15 = 32768 nie powinno być wystarczająco duże, aby wydajność była dużym zmartwieniem.
onestop
Funkcja zwrócona przez splinefun ma domyślnie argument „pochodna”, który wynosi 0. Ustaw pochodną = 1 dla pierwszej pochodnej.
Cyan
1
Hmm, peakswydaje się być wadliwy: Wywołuje max.colz domyślnym ustawieniem ties.method = "random", które nie tylko losowo zrywa więzi, ale także ustawia względną tolerancję 1e-5 dla deklaracji remisu. Pierwsze jest mylące, drugie zdecydowanie nie jest tym, czego chcesz tutaj. peaks()bierze również strictparametr, który jest słabo udokumentowany i patrząc na kod funkcji nic nie robi. Ach, radości z bibliotek oprogramowania tworzonych przez użytkowników! Być może będziesz w stanie to naprawić, ponieważ mówisz, że nie jesteś nowy w programowaniu
onestop

Odpowiedzi:

14

To, co chcesz zrobić, nazywa się wykrywaniem pików w chemometrii. Można do tego użyć różnych metod. Pokazuję tu tylko bardzo proste podejście.

require(graphics)
#some data
d <- density(faithful$eruptions, bw = "sj")

#make it a time series
ts_y<-ts(d$y)

#calculate turning points (extrema)
require(pastecs)
tp<-turnpoints(ts_y)
#plot
plot(d)
points(d$x[tp$tppos],d$y[tp$tppos],col="red")
Roland
źródło
Ze wszystkich rozwiązań działało to najlepiej. 1. Dalsze pytanie: czy istnieje sposób na zmianę tolerancji z punktami zwrotnymi? Znaleziono wiele szczytów i dolin w części ogonowej funkcji zagęszczenia. 2. Dalsze pytanie nr 2: jaki jest dobry sposób ustalenia tolerancji?
aaronlevin
ad 1. Nie sądzę. Jest przeznaczony do testowania losowości szeregów czasowych, więc funkcja tego nie potrzebuje. Możesz spróbować samodzielnie przetestować istotność / znaczenie piku. Np. Możesz wykonać test t względem sąsiedztwa (gdzie możesz zdecydować, jak duże powinno być sąsiedztwo). Możesz też poszukać bardziej wyrafinowanej funkcji w pakietach R do oceny danych ze spektrometrii (masowej) lub innych metod chemii analitycznej.
Roland