Area di bawah "pdf" dalam estimasi kepadatan kernel di R

15

Saya mencoba menggunakan fungsi ' density ' di R untuk melakukan estimasi kepadatan kernel. Saya mengalami beberapa kesulitan menafsirkan hasil dan membandingkan berbagai dataset karena tampaknya area di bawah kurva belum tentu 1. Untuk setiap fungsi kepadatan probabilitas (pdf) , kita perlu memiliki area - ϕ ( x ) d x = 1 . Saya berasumsi bahwa estimasi kepadatan kernel melaporkan pdf. Saya menggunakan integrate.xy dari sfsmisc untuk memperkirakan daerah di bawah kurva.ϕ(x)-ϕ(x)dx=1

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

plot kepadatan

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

kepadatan dengan 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)

kepadatan dengan bw = 1

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

kepadatan dengan bw = 1e-6

Bukankah seharusnya area di bawah kurva selalu 1? Tampaknya bandwidth kecil adalah masalah, tetapi kadang-kadang Anda ingin menunjukkan detail dll di bagian ekor dan bandwidth kecil diperlukan.

Perbarui / Jawab:

220

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

kepadatan dengan jumlah titik sampel yang lebih tinggi pada

> 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

highBandWidth
sumber
3
Ini terlihat seperti batasan titik mengambang dalam kepadatan (): dalam menggunakan bandwidth 1e-6, Anda membuat (secara teori) koleksi 10.000 paku, masing-masing massa total 1/10000. Paku-paku itu akhirnya diwakili terutama oleh puncaknya, tanpa celah yang cukup dikarakterisasi. Anda hanya mendorong kepadatan () di luar batasnya.
whuber
@whuber, dengan batasan floating point, maksud Anda batas presisi, karena dalam menggunakan floats akan menyebabkan perkiraan kesalahan yang lebih besar dibandingkan dengan menggunakan doubles. Saya rasa saya tidak melihat bagaimana itu akan terjadi tetapi saya ingin melihat beberapa bukti.
highBandWidth
n
1
@ Anony-Mousse, ya, itulah yang ditanyakan pertanyaan ini. Mengapa tidak mengevaluasi ke 1?
highBandWidth

Jawaban:

9

Pikirkan tentang penggunaan aturan trapesium integrate.xy(). Untuk distribusi normal, itu akan meremehkan area di bawah kurva dalam interval (-1,1) di mana kepadatannya cekung (dan karenanya interpolasi linier di bawah kepadatan sebenarnya), dan melebih - lebihkannya di tempat lain (saat interpolasi linier berjalan di atas kepadatan sebenarnya). Karena wilayah yang terakhir lebih besar (dalam ukuran Lesbegue, jika Anda suka), aturan trapesium cenderung melebih-lebihkan integral. Sekarang, saat Anda pindah ke bandwidth yang lebih kecil, hampir semua perkiraan Anda sedikit demi sedikit cembung, dengan banyak lonjakan sempit yang terkait dengan titik data, dan lembah di antaranya. Di situlah aturan trapesium rusak parah.

Tugas
sumber
itu berarti bahwa kita "melampauinya" puncak dan "undersampling" lembah-lembah, dalam arti tangan bergelombang. Karena visualisasi juga mengikuti aturan trapesium (interpolasi linier di seluruh sampel), tampaknya terlalu kecil bandwidth kernel juga buruk untuk visualisasi. Juga, jika kita bisa mendapatkan jumlah poin yang lebih besar di mana kita menghitung kepadatan, akan ada sedikit masalah.
highBandWidth
1
Penjelasan ini tidak menahan air. Masalahnya adalah kepadatannya tidak cukup didiskritisasi, bukan karena aturan trapesium rusak parah. integrasikan () tidak berdaya untuk mendapatkan jawaban yang benar karena kepadatan () tidak menghasilkan representasi yang benar. Untuk melihat ini, cukup periksa xy $ x: hanya memiliki nilai 512 yang dimaksudkan untuk mewakili 10.000 paku sempit!
whuber
@whuber, itulah jawabannya. Intinya adalah Anda perlu menggunakan aturan trapesium untuk jumlah sampel yang terbatas, dan itu melebih-lebihkan area dibandingkan dengan kepadatan sebenarnya pada sumbu kontinu sesuai dengan kernel. Pembaruan saya di akhir pertanyaan diperluas di atasnya.
highBandWidth
1
@ Tidak tinggi; aturan trapesium bekerja dengan baik. Masalahnya adalah bahwa ia bekerja dengan diskritisasi yang salah dari integrand. Anda tidak mungkin memiliki "banyak paku sempit yang terkait dengan titik data" ketika ada 10.000 titik data dan hanya 512 nilai dalam array kepadatan!
whuber
1
Melihat grafik-grafik ini, saya sekarang berpikir bahwa masalahnya ada pada densitybukan dengan integrate.xy. Dengan N = 10.000 dan bw = 1e-6, Anda harus melihat sisir dengan tinggi masing-masing gigi sekitar 1e6, dan gigi lebih padat sekitar 0. Sebagai gantinya, Anda masih melihat kurva berbentuk lonceng yang dapat dikenali. Jadi densitymenipu Anda, atau setidaknya itu harus digunakan secara berbeda dengan bandwidth kecil: nharus tentang (rentang data) / (bw) daripada default n=512. Intergrator harus mengambil salah satu dari nilai-nilai besar yang densitydihasilkan oleh suatu kebetulan yang tidak bahagia.
Tugas
-1

Tidak apa-apa, Anda bisa memperbaikinya dengan menggeser dan mengubah skala; tambahkan angka terkecil sedemikian rupa sehingga densitasnya adalah non-negatif, lalu gandakan semuanya dengan konstanta sehingga area tersebut adalah satu. Ini cara yang mudah.

L.2c[ϕ(x)-c]+

Emre
sumber
2
Perhatikan bahwa pertanyaannya adalah lebih pada mengapa yang densityfungsi tidak menghasilkan kepadatan "tepat" yang terintegrasi ke 1 - ketimbang pada bagaimana untuk memperbaikinya.
Tim