“rbinom()”中可能存在大量试验的错误

问题描述 投票:0回答:3

我一直在编写一些迭代执行二项式绘制的代码(使用

rbinom
),对于某些被调用者参数,我最终可能会得到很大的大小,这会导致 R (3.1.1,官方或自制版本都经过测试 - 所以不太可能与编译器相关)返回意外的
NA
。 例如:

rbinom(1,2^32,0.95)

是我期望的工作,但是返回了

NA
。 但是,使用
size=2^31
prob≤0.5
运行是可行的。

精美手册提到当

size < .Machine$integer.max
为 false 时使用反转,这可能是问题所在吗?

r random
3个回答
5
投票

查看

rbinom
对于如此大的尺寸,执行以下等效操作(在C代码中):

qbinom(runif(n), size, prob, FALSE)

确实:

set.seed(42)
rbinom(1,2^31,0.95)
#[1] 2040095619
set.seed(42)
qbinom(runif(1), 2^31, 0.95, F)
#[1] 2040095619

但是:

set.seed(42)
rbinom(1,2^32,0.95)
#[1] NA
set.seed(42)
qbinom(runif(1), 2^32, 0.95, F)
#[1] 4080199349

正如@BenBolker指出的那样,

rbinom
返回一个整数,如果返回值大于
.Machine$integer.max
,例如,大于我的机器上的
2147483647
,则返回
NA
。相反,
qbinom
返回双精度值。我不知道为什么,而且似乎没有记录。

因此,似乎至少存在未记录的行为,您应该报告它。


3
投票

这是预期的行为,但有两个问题: 1)强制引起的NA应该发出警告 2) 应记录离散随机变量具有整数存储模式的事实。

我已经修复了 1),当我有更多时间时,我将修改文档来修复 2)。


3
投票

我同意(在没有文档说明这是一个问题的情况下)这是一个错误。 一个合理的解决方法是使用正态近似,对于如此大的值,这确实应该非常非常好(并且更快)。 (我本来想让这个简短而简单,但最终有点失控。)

rbinom_safe <- function(n,size,prob,max.size=2^31) {
    maxlen <- max(length(size),length(prob),n)
    prob <- rep(prob,length.out=maxlen)
    size <- rep(size,length.out=maxlen)
    res <- numeric(n)
    bigvals <- size>max.size
    if (nbig <- sum(bigvals>0)) {
        m <- (size*prob)[bigvals]
        sd <- sqrt(size*prob*(1-prob))[bigvals]
        res[bigvals] <- round(rnorm(nbig,mean=m,sd=sd))
    }
    if (nbig<n) {
        res[!bigvals] <- rbinom(n-nbig,size[!bigvals],prob[!bigvals])
    }
    return(res)
}

set.seed(101)
size <- c(1,5,10,2^31,2^32)
rbinom_safe(5,size,prob=0.95)
rbinom_safe(5,3,prob=0.95)
rbinom_safe(5,2^32,prob=0.95)

当均值与 0 或 1(以较接近者为准)相差许多标准差时,正态近似应该可以很好地工作。对于大 N 来说这应该没问题,除非 p 非常极端。 例如: n <- 2^31 p <- 0.95 m <- n*p sd <- sqrt(n*p*(1-p)) set.seed(101) rr <- rbinom_safe(10000,n,prob=p) hist(rr,freq=FALSE,col="gray",breaks=50) curve(dnorm(x,mean=m,sd=sd),col=2,add=TRUE) dd <- round(seq(m-5*sd,m+5*sd,length.out=101)) midpts <- (dd[-1]+dd[-length(dd)])/2 lines(midpts,c(diff(sapply(dd,pbinom,size=n,prob=p))/diff(dd)[1]), col="blue",lty=2)

enter image description here

© www.soinside.com 2019 - 2024. All rights reserved.