Подгонка экспоненциального распределения к частотной таблице

У меня есть следующий набор данных:

intervals <- c("0-10", "10-20", "20-30", "30-40", "40-50", "50-75", "75-100", ">100")
int.mean <- c(5.5, 14.3, 24.9, 35.4, 45.2, 63.1, 86.1, NA)
freq <- c(165, 90, 55, 25, 20, 35, 30, 15)

data <- data.frame(intervals, int.mean, freq)

Я пытаюсь подобрать экспоненциальное распределение к данным, чтобы предсказать вероятность того, что значение превысит 150 с определенной степенью уверенности. Я могу подогнать дистрибутив следующим образом:

library(MASS)
fittedexp <- fitdistr(na.exclude(data$int.mean), "exponential")

Однако это не учитывает частоты, поэтому я не уверен, что делаю это правильно. Затем я планирую использовать функцию оптимизации для создания доверительного интервала для предполагаемой вероятности.


r exponential-distribution
person Aesler    schedule 27.02.2020    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Вы имеете дело с категориальной переменной «интервалы», которая создает дискретное наблюдение счетчиков на основе предполагаемой базовой непрерывной переменной, из которой вы взяли точки останова. Какая-то запутанная ситуация с данными. Технически у вас есть данные с интервальной цензурой. Однако, если вы используете экспоненциальное распределение в качестве предположения, те "средние значения", которые вы вычислили, на самом деле являются средними точками, но нельзя ожидать, что они будут средними значениями экспоненциально распределенной переменной. См. ниже мои пересмотренные комментарии о int.means наблюдения. (Так что теперь я расширим свой первоначальный комментарий, включив в него немного кода R.)

Если мы возьмем конечные точки ваших интервалов в качестве переменной разрывов, а также рассчитаем пропорции в наблюдаемых данных, мы получим:

 brks <- c(0, 10,20,30,40,50,75,100,Inf)
 freq <- c(165, 90, 55, 25, 20, 35, 30, 15)
 prop<- freq/sum(freq)
 prop
#-----
[1] 0.37931034 0.20689655 0.12643678 0.05747126 0.04597701 0.08045977 0.06896552 0.03448276
round(prop,2)
[1] 0.38 0.21 0.13 0.06 0.05 0.08 0.07 0.03

Затем мы можем показать, как может «выглядеть» экспоненциально распределенная переменная с аналогичным средним значением (в отношении пропорций), если ее разделить на эти интервалы:

 table( findInterval( rexp(100, 1/15), brks) )/100

   1    2    3    4    5    6    7 
0.47 0.24 0.12 0.08 0.04 0.04 0.01

Итак, мы могли бы попробовать среднее значение выше 15, скажем, 20?

> table( findInterval( rexp(100, 1/20), brks) )/100

   1    2    3    4    5    6    7    8 
0.37 0.24 0.13 0.09 0.07 0.07 0.02 0.01 
> round(prop,2)
[1] 0.38 0.21 0.13 0.06 0.05 0.08 0.07 0.03

Таким образом, вы можете хорошо подобрать нижний предел наблюдений, но экспоненциально распределенная переменная, кажется, имеет несколько «тоньше» хвост. Поскольку вас интересуют верхние границы данных, вы можете захотеть получить лучшее соответствие на верхних границах, но это будет противоречить вашей цели — статистически обоснованному доверительному интервалу. Вы как бы застряли, поскольку ваши данные не являются должным образом «экспоненциальным» набором наблюдений. (Увеличен размер симуляции до 1000, чтобы уменьшить влияние шума.)

> table( findInterval( rexp(1000, 1/25), brks) )/1000

    1     2     3     4     5     6     7     8 
0.329 0.222 0.141 0.103 0.056 0.094 0.034 0.021 
> round(prop,2)
[1] 0.38 0.21 0.13 0.06 0.05 0.08 0.07 0.03

Подгонка там не выглядит ужасной. Если бы параметр скорости экспоненциального распределения был равен 1/25, то это была бы доля наблюдений, превышающая 150:

 1-pexp(150, 1/25)
#[1] 0.002478752

Возможно полезно: http://jsdajournal.springeropen.com/articles/10.1186/s40488-015-0028-6

Вы также можете попробовать выполнить поиск на сайте CrossValidated.com, где существуют некоторые предыдущие обсуждения.

Изменить. Первоначально я думал, что эти значения int.means являются средними точками границ интервала, но это явно не так, поскольку они кажутся близкими к средним точкам, но имеют значительный джиттер вокруг средних точек. Кроме того, эти значения не согласуются с экспоненциальным распределением, так как в наиболее густонаселенном интервале (0-10) наблюдения должны быть «слева» от средней точки, а не даже слева от средней точки. Вероятно, он должен быть 4,0 или 4,5, но уж точно не выше 5,5. Это предполагает, что в основе этого физического процесса лежит какое-то другое распределение, возможно, какое-то гамма-распределение, которое падает до нуля около нуля, но достигает пика в начале интервала 0-10, а затем имеет более длинный хвост.

person IRTFM    schedule 27.02.2020

arrow_upward
1
arrow_downward

Вы можете расширить данные, используя переменную freq, а затем подогнать распределение

data.expand <- data[rep(seq_len(nrow(data)), times=data$freq), ]
head(data.expand, 3); tail(data.expand, 3)

    intervals int.mean freq             intervals int.mean freq
1        0-10      5.5  165        8.12      >100       NA   15
1.1      0-10      5.5  165        8.13      >100       NA   15
1.2      0-10      5.5  165        8.14      >100       NA   15

library(MASS)
with(subset(data.expand, subset=!is.na(int.mean)),
        fitdistr(int.mean,densfun="exponential")
)    

      rate    
  0.041401745 
 (0.002020198)
person Edward    schedule 27.02.2020