Getting Probability Dichte von Daten

Question

Ich muss einige Daten über Internet-Sitzungen für eine DSL-Leitung analysieren. Ich wollte sehen, wie die Sitzungsdauern verteilt sind. Ich dachte, ein einfacher Weg, dies zu tun, würde darin bestehen, ein Wahrscheinlichkeitsdichtediagramm der Dauer aller Sitzungen zu erstellen.

Ich habe die Daten in R geladen und die density() Funktion verwendet. Also, es war so etwas

plot(density(data$duration), type = "l", col = "blue", main = "Density Plot of Duration", 
    xlab = "duration(h)", ylab = "probability density") 

Ich bin neu in R und diese Art von Analyse. Das war, was ich von Google erfahren habe. Ich bekam eine Verschwörung, aber mir blieben einige Fragen. Ist das die richtige Funktion, um das zu tun, was ich versuche, oder gibt es etwas anderes?

In der Grafik fand ich, dass die Y-Achsenskala von 0 ... 1,5 war. Ich verstehe nicht, wie es 1.5 sein kann, sollte es nicht von 0 ... 1 sein?

Auch möchte ich eine glattere Kurve bekommen. Da der Datensatz sehr groß ist, sind die Zeilen wirklich gezackt. Es wäre schöner, wenn sie geglättet werden, wenn ich das vorstelle. Wie würde ich das machen?

Answer 1

Sie sollten mit dem Bandbreitenparameter (bw) herumspielen, um die Glattheit der Kurve zu ändern. Im Allgemeinen macht R einen guten Job und liefert automatisch eine schöne und glatte Kurve, aber das ist vielleicht nicht der Fall für Ihren spezifischen Datensatz.

Wie für den Anruf Sie verwenden, ja, es ist korrekt, type="l" ist nicht erforderlich, es ist die Standardeinstellung zum Plotten von Dichteobjekten. Die Fläche unter der Kurve (d. H. Das Integral von -Inf bis +Inf Ihrer Dichtefunktion) ist = 1.

Nun, ist eine Dichtekurve die beste Sache in Ihrem Fall? Vielleicht, vielleicht nicht ... es hängt wirklich davon ab, welche Art von Analyse du machen willst. Wahrscheinlich wird die Verwendung von hist ausreichend sein, und vielleicht immer informativer, da Sie bestimmte Bins der Dauer auswählen können (siehe ?hist für weitere Informationen).

Answer 2

Wie nico sagte, sollten Sie hist überprüfen, aber Sie können auch die zwei von ihnen kombinieren. Dann könnte man stattdessen die Dichte mit lines aufrufen. Beispiel:

duration <- rpois(500, 10) # For duration data I assume Poisson distributed 
hist(duration, 
    probability = TRUE, # In stead of frequency 
    breaks = "FD",  # For more breaks than the default 
    col = "darkslategray4", border = "seashell3") 
lines(density(duration - 0.5), # Add the kernel density estimate (-.5 fix for the bins) 
    col = "firebrick2", lwd = 3) 

Sollten Sie so etwas wie:

Beachten Sie, dass die Kerndichteschätzung eine Gaußsche Kernel als Standard annimmt. Aber die Bandbreite ist oft der wichtigste Faktor. Wenn Sie density direkt aufrufen meldet es die Standard geschätzte Bandbreite:

> density(duration) 

Call: 
     density.default(x = duration) 

Data: duration (500 obs.);  Bandwidth 'bw' = 0.7752 

     x     y    
Min. : 0.6745 Min. :1.160e-05 
1st Qu.: 7.0872 1st Qu.:1.038e-03 
Median :13.5000 Median :1.932e-02 
Mean :13.5000 Mean :3.895e-02 
3rd Qu.:19.9128 3rd Qu.:7.521e-02 
Max. :26.3255 Max. :1.164e-01 

Hier 0,7752 ist. Überprüfen Sie es auf Ihre Daten und spielen Sie damit herum, wie von nico vorgeschlagen. Vielleicht möchten Sie sich ?bw.nrd ansehen.

Answer 3

Ich wollte dies als Kommentar zur vorherigen Antwort hinzufügen, aber es ist zu groß. Der scheinbare Skew ist auf die Art zurückzuführen, in der die Werte in einem Histogramm zusammengefasst werden. Es ist oft ein Fehler, Histogramme für diskrete Daten zu verwenden. Siehe unten ...

set.seed(1001) 
tmpf <- function() { 
    duration <- rpois(500, 10) # For duration data I assume Poisson distributed 
    hist(duration, 
     probability = TRUE, # In stead of frequency 
     breaks = "FD",  # For more breaks than the default 
     col = "darkslategray4", border = "seashell3", 
     main="",ann=FALSE,axes=FALSE,xlim=c(0,25),ylim=c(0,0.15)) 
    box() 
    lines(density(duration), # Add the kernel density estimate 
     col = "firebrick2", lwd = 3) 
    par(new=TRUE) 
    plot(table(factor(duration,levels=0:25))/length(duration), 
     xlim=c(0,25),ylim=c(0,0.15),col=4,ann=FALSE,axes=FALSE) 
} 

par(mfrow=c(3,3),mar=rep(0,4)) 
replicate(9,tmpf())