Dobre metody dla wykresów gęstości zmiennych nieujemnych w R?

36
plot(density(rexp(100))

Oczywiście cała gęstość na lewo od zera reprezentuje błąd.

Chciałbym podsumować niektóre dane dla statystycznych i chcę uniknąć pytań o to, dlaczego dane nieujemne mają gęstość na lewo od zera. Wykresy służą do sprawdzania losowości; Chcę pokazać rozkład zmiennych według grup leczenia i kontroli. Rozkłady są często wykładnicze. Histogramy są trudne z różnych powodów.

Szybkie wyszukiwanie w Google daje mi pracę statystyczną na nieujemnych jądrach, np .: to .

Ale czy którykolwiek z nich został zaimplementowany w języku R? Czy spośród wdrożonych metod któreś z nich są „najlepsze” w jakiś sposób dla statystyki opisowej?

EDYCJA: nawet jeśli frompolecenie może rozwiązać mój obecny problem, dobrze byłoby wiedzieć, czy ktoś zaimplementował jądra na podstawie literatury na temat nieujemnej oceny gęstości

użytkownik_ogólny
źródło
3
Nie to, o co prosisz, ale nie zastosowałbym oszacowania gęstości jądra do czegoś, co powinno być wykładnicze, szczególnie w przypadku prezentacji dla odbiorców niestatystycznych. Użyłbym wykresu kwantylowo-kwantylowego i wyjaśniłem, że wykres powinien być prosty, jeśli rozkład byłby wykładniczy.
Nick Cox,
6
plot(density(rexp(100), from=0))?
Stéphane Laurent,
4
Jedną rzeczą, którą czasami robiłem dość skutecznie, jest uzyskanie KDE na logach, a następnie przekształcenie oszacowania gęstości (nie zapominając o Jakubie). Inną możliwością byłoby użycie konfiguracji oszacowania gęstości log-splajn, aby wiedział o granicy.
Glen_b
1
Omówiłem metodę transformacji wspomnianą przez @Glen_b w stata-journal.com/sjpdf.html?articlenum=gr0003 (patrz str. 76-78). Zera można dostosować, używając logu (x + 1) zamiast logu i modyfikując jakobian.
Nick Cox

Odpowiedzi:

21

Jednym z rozwiązań, zapożyczonym z podejść do ważenia krawędzi statystyki przestrzennej, jest obcięcie gęstości po lewej stronie na zero, ale zwiększenie wagi danych, które są najbliższe zeru. Chodzi o to, że każda wartość jest „rozkładana” na jądro o jednostkowej całkowitej powierzchni wyśrodkowanej na x ; każda część jądra, która przelałaby się na terytorium ujemne, jest usuwana, a jądro jest normalizowane do obszaru jednostki.xx

Na przykład z jądrem Gaussa , waga renormalizacji wynosiK.h(y,x)=exp(-12)((y-x)/h)2))/2)π

w(x)=1/0K.(y,x)rey=11-Φx,h(0)

gdzie jest funkcją skumulowanego rozkładu normalnej zmiany średniej x i odchylenia standardowego h . Porównywalne formuły są dostępne dla innych jąder.Φxh

Jest to prostsze - i znacznie szybsze w obliczeniach - niż próba zawężenia pasma w pobliżu . W każdym razie trudno jest dokładnie określić, w jaki sposób należy zmienić przepustowość w pobliżu zera . Niemniej jednak ta metoda jest również ad hoc : nadal będzie pewne odchylenie w pobliżu 0 . Wygląda na to, że działa lepiej niż domyślna ocena gęstości. Oto porównanie przy użyciu obszernego zestawu danych:000

Postać

Niebieski pokazuje domyślną gęstość, podczas gdy czerwony pokazuje gęstość skorygowaną dla krawędzi przy . Prawdziwy podstawowy rozkład jest śledzony jako linia kropkowana w celach informacyjnych.0


Kod R.

densityFunkcja w Rskarżą się, że suma wag nie jest jednością, ponieważ chce całkę wszystkich liczb rzeczywistych jako jedność, podczas gdy takie podejście sprawia, że całkę liczby dodatnie równe jedności. Dla sprawdzenia, ta ostatnia całka jest szacowana jako suma Riemanna.

set.seed(17)
x <- rexp(1000)
#
# Compute a bandwidth.
#
h <- density(x, kernel="gaussian")$bw # $
#
# Compute edge weights.
#
w <- 1 / pnorm(0, mean=x, sd=h, lower.tail=FALSE)
#
# The truncated weighted density is what we want.
#
d <- density(x, bw=h, kernel="gaussian", weights=w / length(x))
d$y[d$x < 0] <- 0
#
# Check: the integral ought to be close to 1:
#
sum(d$y * diff(d$x)[1])
#
# Plot the two density estimates.
#
par(mfrow=c(1,1))
plot(d, type="n", main="Default and truncated densities", xlim=c(-1, 5))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
Whuber
źródło
21

Alternatywą jest podejście Kooperberga i współpracowników, oparte na szacowaniu gęstości za pomocą splajnów w celu przybliżenia logarytmicznej gęstości danych. Pokażę przykład wykorzystujący dane z odpowiedzi @ whuber, który pozwoli na porównanie podejść.

set.seed(17)
x <- rexp(1000)

W tym celu musisz zainstalować pakiet logspline ; zainstaluj, jeśli nie jest:

install.packages("logspline")

Załaduj pakiet i oszacuj gęstość za pomocą logspline()funkcji:

require("logspline")
m <- logspline(x)

Poniżej zakładam, że obiekt dz odpowiedzi @ whuber jest obecny w obszarze roboczym.

plot(d, type="n", main="Default, truncated, and logspline densities", 
     xlim=c(-1, 5), ylim = c(0, 1))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
plot(m, add = TRUE, col = "red", lwd = 3, xlim = c(-0.001, max(x)))
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)

Powstały wykres pokazano poniżej, a gęstość logspline jest pokazana czerwoną linią

Domyślne, obcięte i gęstości logspline

Ponadto obsługę gęstości można określić za pomocą argumentów lboundi ubound. Jeśli chcemy założyć, że gęstość wynosi 0 na lewo od 0, a nieciągłość wynosi 0, możemy użyć lbound = 0w wywołaniu logspline()na przykład

m2 <- logspline(x, lbound = 0)

Uzyskano następujące oszacowanie gęstości (pokazane tutaj z oryginalnym mdopasowaniem logspline, ponieważ poprzedni rysunek był już zajęty).

plot.new()
plot.window(xlim = c(-1, max(x)), ylim = c(0, 1.2))
title(main = "Logspline densities with & without a lower bound",
      ylab = "Density", xlab = "x")
plot(m,  col = "red",  xlim = c(0, max(x)), lwd = 3, add = TRUE)
plot(m2, col = "blue", xlim = c(0, max(x)), lwd = 2, add = TRUE)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)
axis(1)
axis(2)
box()

Powstały wykres pokazano poniżej

Porównanie oszacowań gęstości logspline z dolną granicą podpory i bez niej

W tym przypadku wykorzystując wiedzę o xx=0x

Przywróć Monikę - G. Simpson
źródło
1
+1 Podoba mi się ten pomysł, ale jestem zaskoczony, że działa tak dobrze z przykładowymi danymi (z rozkładu wykładniczego). Czy masz intuicję, dlaczego tak jest? W pewnym sensie radzi sobie świetnie w pobliżu ale brakuje w niej „bryły” wartości bliskich 1 w rzeczywistych danych, więc zastanawiam się, czy nie ma jakiegoś kompromisu między dobrą dokładnością przy niskich wartościach i niską dokładnością (lub równoważnie , duże przepustowości) przy wysokich wartościach. 01
whuber
@whuber Dobre pytanie. Takie podejście spotkałem dopiero niedawno. Podejrzewam, że dobrym pytaniem jest tutaj, ponieważ skrócone i logspline metody są jedynie szacunkami rzeczywistej gęstości, czy różnice w dopasowaniu są istotne statystycznie? Nie jestem jednak do końca pewien, dlaczego tak dobrze radzi sobie na zero. Byłbym wdzięczny, wiedząc, dlaczego.
Przywróć Monikę - G. Simpson
@GavinSimpson, Dzięki za tę miłą odpowiedź. Czy możesz odtworzyć ostatnią fabułę w najnowszej wersji logspline? Dla mnie gęstość zarówno wersji ograniczonej, jak i nieograniczonej wynosi zero x = 0.
cel
4

Aby porównać rozkłady według grup (które według ciebie są celem jednego z twoich komentarzy), dlaczego nie coś prostszego? Równoległe wykresy pudełkowe działają dobrze, jeśli N jest duże; wykresy pasków równoległych działają, jeśli N jest małe (i oba pokazują dobrze wartości odstające, co, jak mówisz, jest problemem w twoich danych).

Peter Flom - Przywróć Monikę
źródło
1
Tak, dzięki, to działa. Ale lubię wykresy gęstości. Pokazują więcej o danych niż wykresy. Myślę, że jestem trochę zaskoczony, że wydaje się, że nic już nie zostało wdrożone. Może któregoś dnia zrealizuję jedną z tych rzeczy. Ludzie prawdopodobnie uznaliby to za przydatne.
generic_user
1
Lubię też wykresy gęstości; ale musisz wziąć pod uwagę swoich odbiorców.
Peter Flom - Przywróć Monikę
1
Muszę się zgodzić z @PeterFlom w tej sprawie. Nie komplikuj się zbytnio, jeśli twoi odbiorcy nie posiadają wiedzy statystycznej. Można również wykonać porównawcze / równoległe wykresy pudełkowe z nakładką wykresów motylkowych na górze. W ten sposób podsumowanie wykresu pudełkowego jest widoczne, podobnie jak wszystkie dane.
doug.numbers
Sugestia, że ​​różni ludzie inaczej rozumieją zagregowane wykresy, jest z pewnością słuszna. Pomimo zrozumienia, czym jest wykres gęstości (i zrozumienia, że ​​nie jest to prawdopodobne), nie rozumiem, czym może być „równoległy wykres pudełkowy”. Sugeruje to równoległy wykres współrzędnych, ale podejrzewam, że jest to nieprawidłowe.
DW
2

Jak komentuje Stéphane, możesz użyć, from = 0a dodatkowo możesz przedstawić swoje wartości pod krzywą gęstości za pomocąrug (x)

Aghila
źródło
4
Popraw mnie, jeśli się mylę, ale from=0wygląda na to, że po prostu powstrzymuje wykreślanie wartości poniżej 0; to nie poprawia obliczeń, ponieważ część rozkładu została rozmazana poniżej zera
Nick Cox
1
To jest poprawne. Użycie frompolecenia daje wykres, który wygląda tak, jakby miał szczyt dokładnie na zero. Ale jeśli spojrzysz na histogramy z ciągle mniejszymi pojemnikami, wiele danych pokaże szczytowe AT zero. To fromtylko sztuczka graficzna.
generic_user
@NickCox Nie jestem pewien, ale nie sądzę, aby from=0coś tłumiło. Po prostu zaczyna „siatkę” od zera.
Stéphane Laurent,
Różnica polega na tym, czy szacowana gęstość jest różna od zera dla wartości ujemnych, a nie od tego, czy jest wykreślona. Badacze mogą postanowić nie przejmować się tym, jeśli wszystko, czego chcą, to wizualizacja.
Nick Cox
@NickCox Polecenie density(rexp(100), from=0)nie ma nic wspólnego z grafiką
Stéphane Laurent,