Possible error in `rbinom ()` for a large number of samples

Question

I am writing some code that iteratively executes binomial draws (using rbinom), and for some arguments of the called I can end up having a large size that calls R (3.1.1, both official and homebuilding assemblies - is unlikely to be related with compiler) to return unexpected NA. For instance:

rbinom(1,2^32,0.95)

- this is what I expect to work, but brings NAback. However, it works with size=2^31or prob≤0.5.

The thin guide mentions the inversion used when it size < .Machine$integer.maxis false, could this be a problem?

Answer 1

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 double. , , , .

, , , .

Answer 2

, ( , ), . , ( ) . ( , , - .)

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)![enter image description here][1]
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)

Answer 3

, : 1) , , 2) , .

1) 2), .