Fläche unter dem "pdf" in der Kerndichteschätzung in R

Ich versuche, die ' Dichte' -Funktion in R zu verwenden, um Kernel-Dichteschätzungen durchzuführen. Ich habe einige Schwierigkeiten, die Ergebnisse zu interpretieren und verschiedene Datensätze zu vergleichen, da die Fläche unter der Kurve nicht unbedingt 1 zu sein scheint. Für jede Wahrscheinlichkeitsdichtefunktion (pdf) müssen wir die Fläche . Ich gehe davon aus, dass die Schätzung der Kerneldichte das PDF ausgibt. Ich verwende integrate.xy von sfsmisc die Fläche unter der Kurve zu schätzen. $\phi(x)$ $\int_{-\infty}^\infty \phi(x) dx = 1$

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

Diagramm der Dichte

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

Dichte mit bw = .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)

Dichte mit bw = 1

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

Dichte mit bw = 1e-6

Sollte die Fläche unter der Kurve nicht immer 1 sein? Es scheint, dass kleine Bandbreiten ein Problem sind, aber manchmal möchten Sie die Details usw. in den Endstücken anzeigen, und kleine Bandbreiten werden benötigt.

Update / Antwort:

Es scheint, dass die Antwort über die Überschätzung in konvexen Bereichen richtig ist, da die Erhöhung der Anzahl der Integrationspunkte das Problem zu verringern scheint (ich habe nicht versucht, mehr als Punkte zu verwenden.) $2^{20}$

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

Dichte mit einer höheren Anzahl von Punkten, bei denen eine Probe genommen werden soll

> 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

r estimation pdf kernel-smoothing auc grosse Bandbreite
quelle

Dies sieht aus wie eine Gleitkomma-Beschränkung der Dichte (): Wenn Sie eine Bandbreite von 1e-6 verwenden, erstellen Sie (theoretisch) eine Sammlung von 10.000 Spikes mit einer Gesamtmasse von jeweils 1/10000. Diese Spitzen werden in erster Linie durch ihre Spitzen dargestellt, ohne dass die Lücken ausreichend charakterisiert sind. Sie drücken die Dichte () lediglich über ihre Grenzen hinaus.

Whuber

@whuber, mit Gleitkommabeschränkung meinen Sie Grenzen der Genauigkeit, da bei der Verwendung von Floats der Fehler stärker überschätzt würde als bei der Verwendung von Doubles. Ich glaube nicht, dass ich sehe, wie das passieren würde, aber ich würde gerne Beweise dafür sehen.

HighBandWidth

Ihr Update zeigt, dass Konvexität nicht das Problem ist. Das Problem liegt in der Verwendung eines zu kleinen Wertes von in der Dichteberechnung.

n

$n$

whuber

1

$1$

@ Anony-Mousse, ja, das ist es, was diese Frage stellt. Warum wird nicht 1 ausgewertet?

HighBandWidth

Antworten:

Denken Sie über die integrate.xy()Verwendung der Trapezregel nach . Für die Normalverteilung, werden sie unterschätzen die Fläche unter der Kurve im Intervall (-1,1) , wo die Dichte konkav ist (und somit die lineare Interpolation unter der wahren Dichte) und überschätzen es an anderer Stelle (wie die lineare Interpolation geht auf die wahre Dichte). Da der letztere Bereich größer ist (in Lesbegue, wenn Sie möchten), neigt die Trapezregel dazu, das Integral zu überschätzen. Wenn Sie jetzt auf kleinere Bandbreiten umsteigen, ist fast Ihre gesamte Schätzung stückweise konvex, mit vielen schmalen Spitzen, die den Datenpunkten und Tälern zwischen ihnen entsprechen. Hier bricht die Trapezregel besonders stark zusammen.

StasK
quelle

das bedeutet, dass wir die Gipfel "überabtasten" und die Täler "unterabtasten", in gewissem handgewelltem Sinne. Da die Visualisierung auch der Trapezregel (lineare Interpolation über Samples) folgt, scheint eine zu kleine Kernelbandbreite auch für die Visualisierung schlecht zu sein. Wenn wir eine größere Anzahl von Punkten erhalten könnten, an denen wir die Dichte berechnen, wäre dies weniger problematisch.

HighBandWidth

Diese Erklärung enthält kein Wasser. Das Problem ist, dass die Dichte unzureichend diskretisiert ist und nicht, dass die Trapezregel schlecht zusammenbricht. integriere () ist hilflos, um eine korrekte Antwort zu erhalten, da Dichte () keine korrekte Darstellung erzeugt. Schauen Sie sich dazu einfach xy $ x an: Es gibt nur 512 Werte für 10.000 schmale Spitzen!

whuber

@whuber, das war die Antwort. Der Punkt ist, dass Sie die Trapezregel für eine endliche Anzahl von Samples verwenden müssen und sie die Fläche im Vergleich zur wahren Dichte auf einer kontinuierlichen Achse gemäß den Kerneln überschätzt. Mein Update am Ende der Frage erweitert es.

HighBandWidth

@high Nein; Die Trapezregel funktioniert einwandfrei. Das Problem ist, dass es mit einer falschen Diskretisierung des Integranden arbeitet. Sie können unmöglich "viele schmale Spitzen haben, die den Datenpunkten entsprechen", wenn es 10.000 Datenpunkte und nur 512 Werte im Dichtearray gibt!

Whuber

Wenn ich mir diese Grafiken anschaue, denke ich jetzt, dass das Problem densityeher bei als bei liegt integrate.xy. Mit N = 10000 und bw = 1E6, würden Sie haben einen Kamm mit einer Höhe jedes Zahnes von etwa 1E6 zu sehen, und die Zähne dichter um 0 Stattdessen ist man immer noch ein erkennbares Glockenkurve zu sehen. So densitybetrügt Sie, oder zumindest sollte es anders mit kleinen Bandbreiten verwendet werden: nsollte etwa (Datenbereich) sein / (bw) statt der Standard n=512. Der Intergrator muss einen dieser riesigen Werte densityaufgreifen, die durch einen unglücklichen Zufall zurückkehren.

StasK

-1

Das ist in Ordnung, Sie können es verschieben und skalieren. Addieren Sie die kleinste Zahl, sodass die Dichte nicht negativ ist, und multiplizieren Sie das Ganze mit einer Konstanten, sodass die Fläche eins ist. Das ist der einfache Weg.

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

Emre
quelle

Beachten Sie, dass sich die Frage eher darauf bezieht, warum die densityFunktion nicht die "richtige" Dichte ergibt, die zu 1 integriert wird, als darauf, wie sie zu beheben ist.

Tim