Wydajna obliczeniowo estymacja trybu wielowymiarowego

Wersja skrócona: Jaka jest najbardziej wydajna obliczeniowo metoda szacowania trybu wielowymiarowego zestawu danych, próbkowanego z ciągłego rozkładu?

Wersja długa: Mam zestaw danych, który muszę oszacować dla trybu. Tryb nie pokrywa się ze średnią lub medianą. Przykład pokazano poniżej, jest to przykład 2D, ale rozwiązanie ND byłoby lepsze:

Obecnie moją metodą jest

1. Oblicz oszacowanie gęstości jądra na siatce równej pożądanej rozdzielczości trybu
2. Poszukaj największego obliczonego punktu

Oczywiście oblicza to KDE w wielu niewiarygodnych punktach, co jest szczególnie złe, jeśli istnieje wiele punktów danych o dużych wymiarach lub oczekuję dobrej rozdzielczości w trybie.

Alternatywą byłoby użycie symulowanego wyżarzania, algorytmu genetycznego itp., Aby znaleźć globalny pik w KDE.

Pytanie brzmi, czy istnieje mądrzejsza metoda wykonywania tego obliczenia?

Metodą, która pasuje do rachunku za to, co chcesz zrobić, jest algorytm przesunięcia średniego . Zasadniczo, średni przesunięciem polega na ruchu wzdłuż kierunku gradientu, który ocenia się nie parametrycznie z „cienia”, danego jądra . To znaczy, jeśli gęstość jest szacowana przez , to jest szacowana przez . Szczegóły szacowania gradientu gęstości jądra opisano w tym artykule , który również wprowadził algorytm przesunięcia średniego. K f ( x ) K ∇ f ( x ) K $K'$$K$$f(x)$$K$$\nabla f(x)$$K'$

Bardzo szczegółowy opis algorytmu znajduje się również w tym wpisie na blogu .

Jeśli twoim głównym zainteresowaniem są problemy dwuwymiarowe, powiedziałbym, że oszacowanie gęstości jądra jest dobrym wyborem, ponieważ ma ładne właściwości asymptotyczne (zauważ, że nie twierdzę, że jest najlepszy). Zobacz na przykład

Parzen, E. (1962). O oszacowaniu funkcji i trybu gęstości prawdopodobieństwa . Annals of Mathematical Statistics 33: 1065–1076.

de Valpine, P. (2004). Prawdopodobieństwa przestrzeni stanu Monte Carlo według ważonej oceny gęstości tylnej jądra . Journal of the American Statistics Association 99: 523-536.

W przypadku większych wymiarów (4+) ta metoda jest naprawdę powolna ze względu na dobrze znaną trudność w oszacowaniu optymalnej macierzy przepustowości, patrz .

Problem z poleceniem ksw pakiecie KDEpolega na tym, że, jak wspomniałeś, ocenia on gęstość w określonej siatce, co może być bardzo ograniczające. Ten problem można rozwiązać, jeśli używasz pakietu KDEdo oszacowania macierzy przepustowości, na przykład Hscvzaimplementując estymator gęstości jądra, a następnie optymalizując tę ​​funkcję za pomocą polecenia optim. Jest to pokazane poniżej przy użyciu danych symulowanych i jądra Gaussa R.

rm(list=ls())

# Required packages
library(mvtnorm)
library(ks)

# simulated data
set.seed(1)
dat = rmvnorm(1000,c(0,0),diag(2))

# Bandwidth matrix
H.scv=Hlscv(dat)

# [Implementation of the KDE](http://en.wikipedia.org/wiki/Kernel_density_estimation)
H.eig = eigen(H.scv)
H.sqrt = H.eig$vectors %*% diag(sqrt(H.eig$values)) %*% solve(H.eig$vectors)
H = solve(H.sqrt)
dH = det(H.scv)

Gkde = function(par){
return( -log(mean(dmvnorm(t(H%*%t(par-dat)),rep(0,2),diag(2),log=FALSE)/sqrt(dH))))
}

# Optimisation
Max = optim(c(0,0),Gkde)$par
Max

Na przykład estymatory o ograniczonym kształcie są zwykle szybsze

Cule, ML, Samworth, RJ i Stewart, MI (2010). Oszacowanie maksymalnego prawdopodobieństwa wielowymiarowej gęstości logarytmiczno-wklęsłej . Journal Royal Statistics Society B 72: 545–600.

Ale są one zbyt spiczasty do tego celu.

Problem w dużych wymiarach jest trudny do ataku niezależnie od zastosowanej metody ze względu na charakter samego pytania. Na przykład metoda zaproponowana w innej odpowiedzi (przesunięcie średnie) jest dobra, ale wiadomo, że oszacowanie pochodnej gęstości jest jeszcze trudniejsze niż oszacowanie samej gęstości pod względem błędów (nie krytykuję tego, tylko wskazuję jak trudny jest ten problem). Wtedy prawdopodobnie będziesz potrzebować tysięcy obserwacji, aby dokładnie oszacować tryb w wymiarach większych niż w przypadku problemów innych niż zabawki.$4$

Inne metody, które możesz rozważyć, to: dopasowanie wielowymiarowej skończonej mieszanki normalnych (lub innych elastycznych rozkładów) lub

Abraham, C., Biau, G. i Cadre, B. (2003). Proste oszacowanie trybu wielowymiarowej gęstości . Canadian Journal of Statistics 31: 23–34.

Mam nadzieję, że to pomoże.

Niedawno opublikowaliśmy artykuł sugerujący szybki estymator trybu spójnego.

PS Ruzankin i AV Logachov (2019). Estymator trybu szybkiego w przestrzeni wielowymiarowej. Statystyka i listy prawdopodobieństwa

Nasz estymator ma złożoność czasową , gdzie jest wymiarowością, a jest liczbą obserwowanych punktów. Chociaż nasza metoda może nie być tak precyzyjna, jak inne już tu wspomniane, wypisujemy kompletne dowody na spójność i silną spójność.$O(dn)$$d$$n$

Sugerowałbym również nowe estymatory trybu minimalnej wariancji z mojego ostatniego artykułu

PS Ruzankin (2020). Klasa estymatorów trybu nieparametrycznego. Komunikacja w statystyce - symulacja i obliczenia

Te estymatory mają złożoność czasową dla punktów w . Zobacz rozdział 2.3 tam. Estymatory mają dokładność podobną do znanych algorytmów.$O(dn^2)$$n$${\mathbb R}^d$