Obszar pod „pdf” w szacowaniu gęstości jądra w R

Próbuję użyć funkcji „ gęstości ” w R do oszacowania gęstości jądra. Mam pewne trudności z interpretacją wyników i porównywaniem różnych zestawów danych, ponieważ wydaje się, że obszar pod krzywą niekoniecznie jest 1. Dla każdej funkcji gęstości prawdopodobieństwa (pdf) musimy mieć obszar . Zakładam, że oszacowanie gęstości jądra zgłasza pdf. Korzystam z integrate.xy z sfsmisc, aby oszacować obszar pod krzywą. $\phi(x)$ $\int_{-\infty}^\infty \phi(x) dx = 1$

> # generate some data
> xx<-rnorm(10000)
> # get density
> xy <- density(xx)
> # plot it
> plot(xy)

wykres gęstości

> # load the library
> library(sfsmisc)
> integrate.xy(xy$x,xy$y)
[1] 1.000978
> # fair enough, area close to 1
> # use another bw
> xy <- density(xx,bw=.001)
> plot(xy)

gęstość przy bw = 0,001

> integrate.xy(xy$x,xy$y)
[1] 6.518703
> xy <- density(xx,bw=1)
> integrate.xy(xy$x,xy$y)
[1] 1.000977
> plot(xy)

gęstość z bw = 1

> xy <- density(xx,bw=1e-6)
> integrate.xy(xy$x,xy$y)
[1] 6507.451
> plot(xy)

gęstość przy bw = 1e-6

Czy obszar pod krzywą nie powinien zawsze wynosić 1? Wygląda na to, że małe przepustowości stanowią problem, ale czasami chcesz pokazać szczegóły itp. W ogonach i potrzebne są małe przepustowości.

Aktualizacja / odpowiedź:

$2^{20}$

> xy <- density(xx,n=2^15,bw=.001)
> plot(xy)

gęstość z większą liczbą punktów do próbkowania przy

> integrate.xy(xy$x,xy$y)
[1] 1.000015
> xy <- density(xx,n=2^20,bw=1e-6)
> integrate.xy(xy$x,xy$y)
[1] 2.812398

— wysoka przepustowość
źródło

Wygląda to na zmiennoprzecinkowe ograniczenie gęstości (): używając szerokości pasma 1e-6, tworzysz (teoretycznie) zbiór 10 000 skoków, każdy o masie całkowitej 1/10000. Te kolce są w końcu reprezentowane głównie przez ich szczyty, bez luk odpowiednio scharakteryzowanych. Po prostu przesuwasz gęstość () poza jej granice.

— whuber

@ whuber, przez ograniczenie liczby zmiennoprzecinkowej, masz na myśli granice precyzji, ponieważ użycie liczb zmiennoprzecinkowych doprowadziłoby do większego przeszacowania błędu w porównaniu z użyciem liczb podwójnych. Nie sądzę, że rozumiem, jak to się stanie, ale chciałbym zobaczyć jakieś dowody.

— highBandWidth,

n

$n$

1

$1$

@ Anony-Mousse, tak, właśnie o to pyta. Dlaczego nie ocenia się na 1?

— highBandWidth

Odpowiedzi:

Pomyśl o zasadzie trapezu integrate.xy() zastosowaniach . W przypadku rozkładu normalnego nie doceni on obszaru pod krzywą w przedziale (-1,1), w którym gęstość jest wklęsła (a zatem interpolacja liniowa jest poniżej gęstości rzeczywistej), i zawyża ją w innym miejscu (w miarę interpolacji liniowej na szczycie prawdziwej gęstości). Ponieważ ten drugi obszar jest większy (w miarę Lesbegue, jeśli chcesz), reguła trapezowa ma tendencję do przeceniania całki. Teraz, gdy przechodzisz do mniejszych przepustowości, prawie wszystkie twoje szacunki są częściowo wypukłe, z wieloma wąskimi skokami odpowiadającymi punktom danych i dolinami między nimi. Właśnie tam reguła trapezu psuje się szczególnie źle.

— StasK
źródło

oznacza to, że „nadpróbkujemy” szczyty i „podpróbkujemy” doliny, w pewnym sensie falującym ręcznie. Ponieważ wizualizacja działa również zgodnie z regułą trapezoidalną (interpolacja liniowa między próbkami), wydaje się, że zbyt mała szerokość pasma jądra jest również niekorzystna dla wizualizacji. Gdybyśmy mogli uzyskać większą liczbę punktów, w których obliczamy gęstość, problem byłby mniejszy.

— highBandWidth,

To wyjaśnienie nie zawiera wody. Problem polega na tym, że gęstość jest nieodpowiednio dyskrecjonowana, a nie, że reguła trapezowa źle się psuje. Integrate () nie jest w stanie uzyskać poprawnej odpowiedzi, ponieważ gęstość () nie daje poprawnej reprezentacji. Aby to zobaczyć, po prostu sprawdź xy $ x: ma tylko 512 wartości, które mają reprezentować 10 000 wąskich skoków!

— whuber

@ whuber, tak powiedziała odpowiedź. Chodzi o to, że musisz użyć reguły trapezoidalnej dla skończonej liczby próbek, i przecenia ona powierzchnię w porównaniu z rzeczywistą gęstością na osi ciągłej zgodnie z jądrem. Moja aktualizacja na końcu pytania rozwija się na nim.

— highBandWidth

@high No; zasada trapezoidalna działa dobrze. Problem polega na tym, że działa z niepoprawną dyskretyzacją integrandu. Nie można mieć „wielu wąskich skoków odpowiadających punktom danych”, gdy w tablicy gęstości znajduje się 10 000 punktów danych i tylko 512 wartości!

— whuber

Patrząc na te wykresy, myślę teraz, że problem dotyczy densityraczej niż problemu integrate.xy. Przy N = 10000 i mc = 1e-6, będziesz musiał zobaczyć grzebień o wysokości każdego zęba około 1e6, a zęby będą gęstsze wokół 0. Zamiast tego nadal widzisz rozpoznawalną krzywą w kształcie dzwonu. Tak samo densityoszukuje cię, a przynajmniej powinien być używany inaczej przy małych przepustowościach: npowinien być o (zakres danych) / (mc) niż domyślny n=512. Intergrator musi wychwycić jedną z tych ogromnych wartości, która densitypowraca przez nieszczęśliwy zbieg okoliczności.

— StasK,

-1

W porządku, możesz to naprawić przesuwając i skalując; dodaj najmniejszą liczbę, tak aby gęstość nie była ujemna, a następnie pomnóż całość przez stałą, tak aby obszar był jednością. To jest prosty sposób.

$L_2$ $c$ $\left[\phi(x)-c\right]^+$

— Emre
źródło

Należy zauważyć, że pytanie jest raczej dlaczegodensity funkcja nie powoduje „właściwego” gęstości, który integruje do 1 - a następnie, w jaki sposób to naprawić.

— Tim