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

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

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.$x$$x$

Na przykład z jądrem Gaussa , waga renormalizacji wynosi$K_h(y,x) = \exp(-\frac{1}{2}((y-x)/h)^2) / \sqrt{2\pi}$

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.$\Phi$$x$$h$

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:$0$$0$$0$

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)

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ą

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

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

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).

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)