Как найти отрицательные биномиальные вероятности в R - PullRequest
Новая
       11

Как найти отрицательные биномиальные вероятности в R

0 голосов
/ 02 мая 2018

Я потратил 2 месяца на размышления, достоин ли этот вопрос StackOverflow, и я пришел к выводу, что это так.

Я добровольно работаю в команде в течение года, чтобы предсказать ряд интересных вещей, несколько месяцев назад она предсказывала вероятность количества землетрясений во всем мире в течение 5 месяцев в марте. Действительно интересная проблема. Думал, что я неплохо справился с R, а затем столкнулся с этой проблемой, как кирпичная стена. Это проблема подсчета, я хотел бы использовать распределение Пуассона, но это не сработает, среднее и дисперсия не равны. Он рассредоточен.

Цель состоит в том, чтобы оценить вероятность:

<100 землетрясений 100-140 землетрясений 140-170 землетрясений 170-210 <br> землетрясения 210 землетрясений

Но я написал здесь немного кода: 

#(load data and libaries blah blah blah)
quakes_this_month<-10
days_left=31-1
days_left
month_left<- days_left/31
month_left
earthq5<- earthq4
earthq5$mag<-earthq5$mag*month_left
mu <- mean(earthq5$mag)
sigma <- sd(earthq5$mag)
paste("mean is ",mu, " and sigma is ", sigma)
pnorm((99-quakes_this_month) , mu, sigma, lower.tail = T)
lower.bound<- 100 -quakes_this_month
upper.bound<- 140.5-quakes_this_month
(pnorm(upper.bound, mu, sigma, lower.tail = T) - pnorm(lower.bound, mu, sigma))
lower.bound<- 140.5-quakes_this_month
upper.bound<- 170.5-quakes_this_month
(pnorm(upper.bound, mu, sigma) - pnorm(lower.bound, mu, sigma))
lower.bound<- 170.5-quakes_this_month
upper.bound<- 210.5-quakes_this_month
(pnorm(upper.bound, mu, sigma) - pnorm(lower.bound, mu, sigma))
(pnorm(210.5-quakes_this_month, mu, sigma, lower.tail = F))

Так что идея здесь в том, что месяц прогрессирует, и произошло несколько землетрясений, и я могу оценить вероятность достижения этих предельных порогов. Однако это не гауссовский дистрибутив, я не могу использовать pnorm, я должен использовать pnbinom(q, size, prob, mu, lower.tail = TRUE, log.p = FALSE), но я не знаю, как вывести «size» и «prob» из проблемы с подсчетом. Это не берет 15 шаров из банки 4 раза. Так что я обращаюсь к этому, так как он преследовал меня в течение нескольких недель. Как я могу использовать pnbinom() вместо pnorm(), если речь идет о количестве землетрясений в месяц?

1 Ответ

0 голосов
/ 07 октября 2018

Итак, я нашел ответ, и для всех остальных вот как я это сделал. Данные, которые я использовал, были из Геологической службы США о землетрясениях. В R. Я использую довольно много других библиотек. Я думаю, что для этого примера требуется только MASS.

Загрузка библиотеки и данных 

library(MASS)

quakeSim <-  earthq4$count  # this was my real data

quakeSim <-  rnbinom(n = 12000, mu = 145, size =18)  # you can use this for the example

Тест на соответствие распределения распределения 3 вероятных распределения: гауссовское, пуассоновское и отрицательное биномиальное 

  quakeDistNB <- MASS::fitdistr(quakeSim, densfun = "negative binomial")
    quakeDistPois <- MASS::fitdistr(quakeSim, densfun = "poisson")
    quakeDistGaus<-MASS::fitdistr(quakeSim, densfun = "normal")

Сравните отрицательный бином, Пуассона и Гуассиана - чем ниже AIC, тем лучше выберите распределение с самым низким AIC. 

 AIC(quakeDistNB)
    AIC(quakeDistPois)
    AIC(quakeDistGaus)

Быстрая проверка Normalicy с помощью теста Шапиро. (если гауссовский самый низкий) 

shapiro.test(earthq4$count)

Используйте правило 5%. Но это NB, а не гауссов, так что игнорируйте все тесты CI ниже 

summary(earthq4)
t.test(earthq4$count ) #default 0.95

Итак, мои данные показывают отрицательное биномиальное распределение. Теперь давайте посмотрим на это как на гистограмму с достаточным количеством столбцов, чтобы показать форму NB.

визуализация эмпирического распределения 

hist(quakeSim, breaks=80)

Установите отрицательную биномиальную модель и получите два критических значения sizeHat и muHat из выходных данных модели 'quakeDistNB'

Эта часть действительно сводила меня с ума, пока друг не показал мне. 

quakeDistNB <- MASS::fitdistr(earthq4$count , densfun = "negative binomial")
quakeDistNB
sizeHat <- quakeDistNB$estimate[1]
sizeHat
muHat <- quakeDistNB$estimate[2]

Теперь моя задача состояла в том, чтобы предсказать вероятность менее 100 землетрясений и от 150 до 100 с магнитудой больше или равной 5.

Тогда вероятность меньше 100: 

p100 <- pnbinom(q = 100, size = sizeHat, mu = muHat)
p100

вероятность менее 150: 

p150 <- pnbinom(q = 150, size = sizeHat, mu = muHat)
p150

вероятность от 100 до 150: 

p150 - p100
...