Подгонка экспоненциального распределения к таблице частот - PullRequest
Новая
       32

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

0 голосов
/ 27 февраля 2020

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

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

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

Ответы [ 2 ]

1 голос
/ 27 февраля 2020

Вы имеете дело с категориальной переменной, «интервалами», которая создает дискретное наблюдение за счетами, основанными на предполагаемой базовой непрерывной переменной, из которой вы взяли контрольные точки. Вроде запутанная ситуация с данными. Технически у вас есть данные с интервальной цензурой. Однако, если вы используете экспоненциальное распределение в качестве допущения, то рассчитанные вами «средние» на самом деле являются средними точками, но не следует ожидать, что они будут средством экспоненциально распределенной переменной. См. Мои пересмотренные комментарии ниже о 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, а затем имело более длинный хвост.

1 голос
/ 27 февраля 2020

Вы можете расширить данные, используя переменную 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)
...