Bug 17155

Summary: rbeta function returns fixed values when any shape parameter is NA
Product: R Reporter: Joan Maspons <joanmaspons>
Component: MiscAssignee: R-core <R-core>
Status: CLOSED FIXED    
Severity: minor CC: joanmaspons, maechler
Priority: P5    
Version: R 3.3.*   
Hardware: All   
OS: Linux   
Attachments: Patch

Description Joan Maspons 2016-09-20 09:53:23 UTC
Created attachment 2159 [details]

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)
Comment 1 Joan Maspons 2016-09-20 10:20:39 UTC
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="_"),
Comment 2 Joan Maspons 2016-09-20 15:18:20 UTC
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),

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))
Comment 3 Martin Maechler 2016-09-30 07:40:05 UTC
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