Created attachment 2159 [details] Patch rbeta function returns a fixed numeric value when any shape parameter is NA: rbeta(n=5, shape1=1, shape2=1) # [1] 0.3810525 0.3450395 0.8499018 0.4855583 0.9808052 rbeta(n=5, shape1=1, shape2=NA) # [1] 0 0 0 0 0 0 0 0 0 0 rbeta(n=5, shape1=NA, 1) # [1] 1 1 1 1 1 1 1 1 1 1 rbeta(n=5, shape1=NA, shape2=NA) The documentation says: «Invalid arguments will result in return value NaN, with a warning.» But in fact d/p/qbeta functions return NA_real_ instead of NaN. I would expect that rbeta returns a vector of length n filled with NA_real_. This is the same behavior as d/p/qbeta functions. Patch attached The documentation should be updated (NaN -> NA)

A test for the patch and beta distribution functions in general: p<- c(1, NA, Inf) p<- expand.grid(shape1=shape, shape2=shape) dbeta(.5, shape1=p$shape1, shape2=p$shape2) pbeta(.5, shape1=p$shape1, shape2=p$shape2) qbeta(.5, shape1=p$shape1, shape2=p$shape2) replicates<- 10 matrix(rbeta(nrow(p) * replicates, shape1=p$shape1, shape2=p$shape2), nrow=nrow(p), ncol=replicates, dimnames=list(shape1_shape2=paste(p$shape1, p$shape2, sep="_"), replicate=NULL))

Some error in the test. It should be p<- c(1, NA, Inf) p<- expand.grid(shape1=p, shape2=p) pars<- paste(p$shape1, p$shape2, sep="_") t(data.frame(dbeta=dbeta(.5, shape1=p$shape1, shape2=p$shape2), pbeta=pbeta(.5, shape1=p$shape1, shape2=p$shape2), qbeta=qbeta(.5, shape1=p$shape1, shape2=p$shape2), row.names=pars)) replicates<- 10 matrix(rbeta(nrow(p) * replicates, shape1=p$shape1, shape2=p$shape2), byrow=TRUE, ncol=nrow(p), nrow=replicates, dimnames=list(replicate=NULL, shape1_shape2=pars))

Thank you, Joan. Indeed, these boundary cases should be fixed. To get consistency here, I'd like to look at other distributions, not *beta alone, e.g. the normal > rnorm(4, NA) [1] NaN NaN NaN NaN Warning message: In rnorm(4, NA) : NAs produced > This corresponds to what the documentation says for rbeta(), too, and I would want to adapt the behavior to the documentation (which most probably once was correct). Note that this is also what runif(), rt(), and rchisq() produce in analogous examples. We definitely want rbeta() to behave the same here. --- BTW: I do note that rgamma() has a similar issue and I assume there will more such cases.... though they are not "independent". ---> I have written testing-code and fixed all three cases I found: rbeta(), rgamma() and rnbinom(). ---> Fixed in R-patched and R-devel