Получение вероятности плотности данных

Question

Мне нужно проанализировать некоторые данные об интернет-сессиях для линии DSL.Я хотел посмотреть, как распределяются длительности сеанса.Я подумал, что простой способ сделать это - начать с построения графика плотности вероятности продолжительности всех сеансов.

Я загрузил данные в R и использовал функцию density().Итак, это было что-то вроде

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

Я новичок в R и этот вид анализа.Это было то, что я нашел, пройдя через Google.Я получил сюжет, но у меня остались некоторые вопросы.Это правильная функция, чтобы делать то, что я пытаюсь сделать, или есть что-то еще?

На графике я обнаружил, что масштаб оси Y был от 0 ... 1,5.Я не понимаю, как это может быть 1,5, не должно ли это быть от 0 ... 1?

Кроме того, я хотел бы получить более плавную кривую.Так как набор данных действительно большой, строки действительно зазубренные.Было бы лучше, чтобы они были сглажены, когда я представляю это.Как мне поступить так?

Answer 1

Как сказал Нико, вы должны проверить hist, но вы также можете объединить их. Тогда вы могли бы назвать плотность с lines вместо. Пример:

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)

Должно дать вам что-то вроде:

Обратите внимание, что для оценки плотности ядра по умолчанию используется ядро ​​Гаусса. Но пропускная способность часто является наиболее важным фактором. Если вы позвоните density напрямую, он сообщит о предполагаемой пропускной способности по умолчанию:

> 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  

Здесь это 0,7752. Проверьте это на наличие данных и поиграйтесь с ним, как предложил Нико. Возможно, вы захотите посмотреть на ?bw.nrd.

Answer 2

Вам следует поиграться с параметром bandwith (bw), чтобы изменить плавность кривой. Обычно R хорошо выполняет свою работу и автоматически дает хорошую и плавную кривую, но, возможно, это не относится к вашему конкретному набору данных.

Что касается вызова, который вы используете, да, это правильно, type="l" не требуется, это значение по умолчанию, используемое для построения объектов плотности. Площадь под кривой (то есть интеграл от -Inf до + Inf вашей функции плотности) будет = 1.

Теперь, является ли кривая плотности лучшим вариантом для использования в вашем случае? Может быть, может и нет ... это действительно зависит от того, какой анализ вы хотите сделать. Возможно, будет достаточно использовать hist, а может быть, даже более информативно, поскольку вы можете выбрать конкретные интервалы длительности (см. ?hist для получения дополнительной информации).

Answer 3

Я собирался добавить это как комментарий к предыдущему ответу, но он слишком большой. Очевидный перекос связан с тем, как значения объединены в гистограмме. Часто ошибочно использовать гистограммы для дискретных данных. Смотри ниже ...

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