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

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?

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.

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 .

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.

Możesz spróbować zmniejszyć przepustowość (niebieska linia oznacza adjust=0.5),

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