Странные проблемы с точностью в R при вычислении кумулятивной биномиальной вероятности

Я столкнулся с некоторыми странными проблемами при использовании этого кода:

positions<-c(58256)
occurrencies<-c(30)
frequency<-c(11/5531777)
length<-c(4)

prob<-c(0)
for(i in 0:(occurrencies-1))
{
  pow<-frequency^i
  pow1<-(1-frequency)^(positions-i)
  bin<-choose(positions, i)
  prob<<-prob+(bin*pow*pow1)
}

Каждая итерация этого цикла for должна вычислять биномиальную вероятность того, что i количество вхождений события произойдет с заданной частотой. Каждая итерация также подводит итог. Это должно привести к тому, что переменная prob никогда не превысит 1, но после 7 или около того итераций цикла все пойдет к черту, и prob превысит 1.

Я подумал, что это может быть вопрос точности цифр, поэтому я попытался использовать Rmpfr, но безрезультатно - та же проблема осталась.

Мне было интересно, есть ли какие-либо советы или пакеты для преодоления этой ситуации, или я застрял с этим.


r binomial-coefficients arbitrary-precision
person Tito Candelli    schedule 11.10.2012    source источник
comment
Я скопировал/вставил ваш код, и хотя prob действительно больше 1, это ненамного. abs(prob - 1) < 1e-10 возвращает TRUE. Я могу получить дополнительную цифру, изменив последнюю строку на prob <- prob + exp(log(bin) + log(pow) + log(pow1)) (обратите внимание, что <<- не требуется, если вы не заключаете это в функцию и не хотите return(prob)), но если вы хотите большего, просто используйте pbinom как предлагает Бен.   -  person Gregor Thomas    schedule 11.10.2012
comment
Побочный вопрос: когда вы использовали Rmpfr , какой уровень точности вы выбрали, и уверены ли вы, что вычислили frequency как отношение class:mpfr чисел?   -  person Carl Witthoft    schedule 11.10.2012
comment
@shujaa: это правда, что он ненамного больше 1, но меня беспокоит то, что после определенного количества циклов цикла for prob остается прежним, даже если что-то добавляется: попробуйте запустить цикл for от 0 до вхождения-3, чтобы проверить, что это правда. что касается побочного вопроса: я этого не делал, частота на самом деле является числом Rmpfr, происходящим из соотношения двух нормальных чисел. я постараюсь держать вас в курсе.   -  person Tito Candelli    schedule 11.10.2012

Ответы (2)


arrow_upward
2
arrow_downward

Следуя совету Бена Болкера, посмотрите ?pbinom

pbinom(q = occurencies, size = positions, prob = frequency, lower.tail = FALSE)
person Gregor Thomas    schedule 12.10.2012

arrow_upward
3
arrow_downward

Вы можете избежать цикла for, выполнив

prob<-0
i    <- 0:(occurrencies-1)
pow  <- frequency^i
pow1 <- (1-frequency)^(positions-i)
bin  <- choose(positions, i)
prob <- cumsum(prob+(bin*pow*pow1))
[1] 0.8906152 0.9937867 0.9997624 0.9999932 0.9999998 1.0000000 1.0000000 1.0000000 1.0000000
[10] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[19] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[28] 1.0000000 1.0000000 1.0000000

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

См. комментарий @Ben Bolker и взгляните на функцию pbinom.

person Jilber Urbina    schedule 11.10.2012
comment
это именно то, чего я пытаюсь избежать. я не знал о функции cumsum, и я приветствую подсказку; тем не менее, я пытаюсь вычислить вероятность того, что с учетом частоты событие происходит меньше, чем раз, так что это никогда не должно быть 1. что мне действительно нужно, так это вероятность того, что событие происходит или больше раз, но вычисление этих чисел просто сумасшедший, поэтому я решил вычислить обратное. - person Tito Candelli; 11.10.2012
comment
@user1738815 user1738815 Похоже, вы ищете pbinom(q = occurencies, size = positions, prob = frequency, lower.tail = FALSE) - person Gregor Thomas; 12.10.2012
comment
Большое спасибо @shujaa! это именно то, что мне было нужно! как выбрать ваш ответ как правильный? я могу только проголосовать - person Tito Candelli; 12.10.2012
comment
Я добавлю это как ответ, тогда это будет приемлемо. Но на самом деле это просто повторение Бена Болкера. Для справки в будущем, мы могли бы помочь вам только тогда, когда вы сказали нам, что мне действительно нужно в комментарии выше. - person Gregor Thomas; 12.10.2012
comment
Я буду иметь это в виду для будущих вопросов. Еще раз спасибо. - person Tito Candelli; 13.10.2012