To jest tylko wersja do druku, aby zobaczyć pełną wersję tematu, kliknij TUTAJ
Statystyka, prognozowanie, ekonometria, data mining
Forum miłośników statystyki - Portal Statystyczny

Giełda ogłoszeń, nawiązanie współpracy, potrzebna pomoc - Kernel Density Estimation dla R

GZabinski - 2018-01-15, 14:31
Temat postu: Kernel Density Estimation dla R
Dzień dobry,
Chciałbym zapytać, czy ktoś z użytkowników mógłby coś zasugerować w sprawie modyfikacji lub napisania od nowa kodu do R, obliczającego KDE dla wielu zmiennych. Oryginalny kod pochodzi stąd:
https://biostatmatt.com/archives/2745

Oblicza on KDE dla dwóch zmiennych i pokazuje wykres. Moje dane (są to wyniki PCA z badań składu chemicznego wtrąceń żużla) są nieco bardziej złożone, mają następującą postać:
F1 F2 F3 F4 F5 F6 Region
-0.371 0.128 -0.130 0.439 -0.081 0.028 DenmarkWest
-0.357 0.148 0.143 0.176 -0.190 0.034 DenmarkWest
i tak dalej, w sumie 118 obserwacji dla 4 regionów

Dane te były już analizowane wcześniej w literaturze za pomocą KDE (niestety, bez podania szczegółów technicznych, nie mówiąc o kodzie) i otrzymano wynik pokazujący rozdzielenie między regionami.
Przeanalizowałem te dane metodą Discriminant Analysis i otrzymuję dość wyraźne rozdzielenie na obserwacje z poszczególnych regionów. Chciałbym to sprawdzić za pomocą KDE w R, stąd zmodyfikowałem powyższy kod jak następuje:

## radially symmetric kernel (Gaussian kernel)
RadSym <- function(u)
exp(-rowSums(u^2)/2) / (2*pi)^(ncol(u)/2)

## multivariate extension of Scott's bandwidth rule
Scott <- function(data)
t(chol(cov(data))) * nrow(data) ^ (-1/(ncol(data)+4))

## compute KDE at x given data
mvkde <- function(x, data, bandwidth=Scott, kernel=RadSym) {
# bandwidth may be a function or matrix
if(is.function(bandwidth))
bandwidth <- bandwidth(data)
u <- t(solve(bandwidth, t(data) - x))
mean(kernel(u))
}

## compute KDE at (matrix) x given data
smvkde <- function(x, ...)
apply(x, 1, mvkde, ...)

# make a dataset “Slag” with the following structure
Slag= read.delim("clipboard",header=TRUE) #copy a table from Excel, all “,” must be replaced by “.”

#introduce division into data and labels

## compute bivariate KDE and plot contours
data("Slag")
aq <- subset(Slag, !is.na(F1) & !is.na(F2),
select=c("F1", "F2"))

## compute density on a grid of F1 and F2 values
## -0.5 and +0.5 added to densities for a better plot
dens.F1 <- seq(min(aq$F1)-0.5,max(aq$F1)+0.5,length.out=100)
dens.F2 <- seq(min(aq$F2)-0.5,max(aq$F2)+0.5,length.out=100)
dens.grid <- expand.grid(F1=dens.F1, F2=dens.F2)
dens.vals <- smvkde(dens.grid, data=aq)

## arrange density values into matrix for easy plotting
dens.mtrx <- matrix(dens.vals, 100, 100)
contour(x=dens.F1, y=dens.F2, z=dens.mtrx,
xlab="F1", ylab="F2")
points(aq$F1, aq$F2, pch=20)

Niestety, wynikiem jest graf, gdzie wszystkie punkty są na jednym obszarze. Stąd moje pytanie - czy ten kod da się jakoś zmodyfikować, czy należy napisać go na nowo? Brak rozdzielenia między grupami wedle mnie wynika z tego, że kod bierze pod uwagę tylko F1 i F2, a ja chciałbym, żeby uwzględnił wszystko (F1-F2) i pokazał wynik w postaci konturów na osiach F1 i F2.
Artykuł zaleca jeszcze inną metodę wygładzania – w kodzie jest Scott, a ma być Sheather-Jones plug in method.

Ponieważ dopiero zaczynam zabawę z R, będę zobowiązany za pomoc.
Z ukłonami
G. Żabiński

Crunchy - 2018-01-16, 13:20

Sprawdź co jest w dens.vals, jeżeli coś dziwnego, to jest to problem numeryczny. Zdecydowanie lepiej jest korzystać z gotowych (czyt. zwykle lepiej przetestowanych) funkcji z pakietów, kde {ks}, sjpi {locfit}.
GZabinski - 2018-01-16, 13:59

Dzięki, będę próbował iść tym torem
GZabinski - 2018-01-20, 00:49

No dobra, zgodnie z radą Crunchy'ego korzystam z funkcji kde w ks. Ma być dla 1 do 6 wymiarów, ale powyżej 4 trzeba podać eval.points jako matrycę. Ktoś wie, jak to zrobić?
Pokłony
GZ

Crunchy - 2018-01-20, 16:42

Wie, tylko nie wie po co chcesz to zrobić, bo niby czytasz dokumentację, ale...
Cytat:
–If eval.points is specified, as a vector (d=1) or as a matrix (d=2, 3, 4), then the density estimate is computed at eval.points. This form is suitable for numerical summaries (e.g. maximum likelilood), and is not compatible with the plot method.
–For d>4, computing the kernel density estimate over a grid is not feasible, and so it is computed exactly and eval.points (as a matrix) must be specified.
Plot for kernel density estimate for 1- to 3-dimensional data.

GZabinski - 2018-01-20, 18:05

To jak wie, to niech będzie tak uprzejmy i powie ;-)
Ten akapicik oczywiście przeczytałem. Rozumiem, że pod Twoim "niby...., ale...." kryje się sprawa plotu, możliwego jedynie dla max 3 wymiarów?
Szczerze mówiąc, nie wiem, w jaki dokładnie sposób ten wykres (jakoś nie mogę wrzucić obrazka z artykułu, gdzie zastosowano tę metodę) zrobiono, ale autorzy twierdzili, że przeanalizowali metodą KDE zbiór z 6 zmiennymi, a potem wrzucili to na wykres 2 wymiarowy z osiami F1-F2. Staram się zatem dość do tego, jaką metodę zastosowali.
Poza tym, jak zaznaczyłem wyżej, nie znam R prawie wcale (zainstalowałem to 3 tygodnie temu). Nie dziw się zatem, że to i owo nie do końca musi być dla mnie jasne.

Crunchy - 2018-01-20, 19:37

Nie ma takiej możliwości, żeby sensownie przedstawić sześć wymiarów na statycznym wykresie dwuwymiarowym. I to by było na tyle, w tej kwestii.
Musisz mieć dwa wymiary. W jaki sposób to osiągniesz, to już inna bajka, bo możliwości jest wiele.

GZabinski - 2018-01-20, 22:05

To inaczej postawię problem, przez analogię do LDA. Analizuję tą metodą mój zbiór danych z 6 zmiennymi, dostaję wyniki F1-F6. Na F1 i F2 razem mam ponad 90%. Tworzę wykres dwuwymiarowy F1-F2 i dostaję wyraźny podział na 4 grupy regionalne. Pytanie, jak podobną operację dokonać za pomocą KDE w R. Różnica między LDA i KDE ma polegać na tym, że LDA rysuje granice zbiorów w postaci elips, a KDE ma to zrobić inaczej, bliżej punktów poszczególnych wyników. Pozwoli to na uzyskanie lepszego rozdzielenia między zbiorami.
Crunchy - 2018-01-21, 10:29

Wystarczy użyć odpowiadających funkcji. To, jakie kontury zaznaczy możesz sobie ustawić wedle uznania. Jeżeli chcesz mieć jak najlepszą klasyfikację, musisz sprawdzić różne klasyfikatory.
GZabinski - 2018-01-27, 23:32

Rozwiązanie zasadniczo jest. Można ewentualnie jeszcze dopracować. Niestety, obrazek dokumentujący efekt nie pozwala się załadować
ukłony
GZ ;-)



Powered by phpBB modified by Przemo © 2003 phpBB Group