Jak mogę oszacować gęstość parametru z napompowaniem zerowym w R?

10

Mam zestaw danych z dużą ilością zer, który wygląda następująco:

set.seed(1)
x <- c(rlnorm(100),rep(0,50))
hist(x,probability=TRUE,breaks = 25)

Chciałbym narysować linię dla jej gęstości, ale density()funkcja wykorzystuje ruchome okno, które oblicza ujemne wartości x.

lines(density(x), col = 'grey')

Istnieją density(... from, to)argumenty, ale wydają się one jedynie przycinać obliczenia, a nie zmieniać okna, tak aby gęstość przy 0 była zgodna z danymi, co widać na poniższym wykresie:

lines(density(x, from = 0), col = 'black')

(gdyby interpolacja została zmieniona, oczekiwałbym, że czarna linia miałaby większą gęstość przy 0 niż szara linia)

Czy istnieją alternatywy dla tej funkcji, które zapewniłyby lepsze obliczenie gęstości przy zeru?

wprowadź opis zdjęcia tutaj

Abe
źródło

Odpowiedzi:

14

Gęstość jest nieskończona przy zera, ponieważ zawiera dyskretny skok. Musisz oszacować wartość szczytową za pomocą proporcji zer, a następnie oszacować dodatnią część gęstości, zakładając, że jest ona gładka. KDE spowoduje problemy na lewym końcu, ponieważ przywróci pewną wartość ujemnym. Jednym z przydatnych podejść jest przekształcenie w dzienniki, oszacowanie gęstości za pomocą KDE, a następnie przekształcenie z powrotem. Zobacz Wand, Marron i Ruppert (JASA 1991) w celach informacyjnych.

Następująca funkcja R wykona transformowaną gęstość:

logdensity <- function (x, bw = "SJ") 
{
    y <- log(x)
    g <- density(y, bw = bw, n = 1001)
    xgrid <- exp(g$x)
    g$y <- c(0, g$y/xgrid)
    g$x <- c(0, xgrid)
    return(g)
}

Następnie następujące czynności podadzą pożądaną fabułę:

set.seed(1)
x <- c(rlnorm(100),rep(0,50))
hist(x,probability=TRUE,breaks = 25)
fit <- logdensity(x[x>0]) # Only take density of positive part
lines(fit$x,fit$y*mean(x>0),col="red") # Scale density by proportion positive
abline(v=0,col="blue") # Add spike at zero.

wprowadź opis zdjęcia tutaj

Rob Hyndman
źródło
P.(X=0)
P.(X=0)
to się przydaje. fyi: wydaje się, że chociaż bw = „SJ” wpływa na gęstość w przestrzeni nietransformowanej, logdęstość jest taka sama przy użyciu „SJ” i domyślnego „nrd0” ... Zaraz przeczytam odniesienie SJ: „Sheather i Jones (1991) Niezawodna metoda wyboru przepustowości oparta na danych do oceny gęstości jądra. ” jstor.org/stable/2345597
Abe
4

Zgadzam się z Robem Hyndmanem, że musisz zajmować się zerami osobno. Istnieje kilka metod radzenia sobie z estymacją gęstości jądra zmiennej z ograniczonym wsparciem, w tym „odbicie”, „ponowna normalizacja” i „kombinacja liniowa”. Nie wydaje się, aby zostały zaimplementowane w densityfunkcji R , ale są dostępne w pakiecie Benna Janna kdensdla Staty .

jeden przystanek
źródło
1

Inna opcja, gdy masz dane z logiczną dolną granicą (takie jak 0, ale mogą to być inne wartości), o których wiesz, że dane nie spadną poniżej, a oszacowanie normalnego zagęszczenia jądra umieszcza wartości poniżej tej granicy (lub jeśli masz górną granicę lub oba) to użycie oszacowań logspline. Pakiet logspline dla R implementuje je, a funkcje mają argumenty do określenia granic, więc oszacowanie przejdzie do granicy, ale nie dalej i nadal będzie skalowane do 1.

Istnieją również metody ( oldlogsplinefunkcja), które wezmą pod uwagę cenzurowanie interwałów, więc jeśli te 0 nie są dokładnie zerami 0, ale są zaokrąglone, abyś wiedział, że reprezentują wartości między 0 a jakąś inną liczbą (na przykład limit wykrywania) może przekazać tę informację funkcji dopasowania.

Jeśli dodatkowe 0 to prawdziwe 0 (nie zaokrąglone), wówczas oszacowanie wartości szczytowej lub masy punktowej jest lepszym podejściem, ale można je również połączyć z estymacją logspline.

Greg Snow
źródło
0

Możesz spróbować zmniejszyć przepustowość (niebieska linia oznacza adjust=0.5), wprowadź opis zdjęcia tutaj

ale prawdopodobnie KDE nie jest najlepszą metodą radzenia sobie z takimi danymi.


źródło
czy jest jakaś inna metoda, którą poleciłbyś?
Abe
@ Cóż, to zależy od tego, co chcesz zrobić ...