Bug 17375 - rbeta() and rf() don't handle NA in ncp properly
Summary: rbeta() and rf() don't handle NA in ncp properly
Status: CLOSED FIXED
Alias: None
Product: R
Classification: Unclassified
Component: Models (show other bugs)
Version: R 3.4.1
Hardware: Other Other
: P5 enhancement
Assignee: R-core
URL:
Depends on:
Blocks:
 
Reported: 2018-01-11 10:21 UTC by Duncan Murdoch
Modified: 2018-01-11 16:06 UTC (History)
1 user (show)

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Duncan Murdoch 2018-01-11 10:21:17 UTC
Received this from Jonas Moss <jonasmgj@math.uio.no>.  I tested it in R-devel and still see the same thing.  I've modified his suggested fix for rbeta slightly.

(Using "R version 3.4.3 (2017-11-30)").

Example of bad behaviour for rbeta:

set.seed(313)
n = 10
shape1 = 5
shape2 = 20
ncp = c(1:9, NA)
rbeta(n, shape1, shape2, ncp) # [1] 0.2714894 0.2753471 0.2715325 0.2736262 0.2349107 0.4290410 0.2594788 0.1888690 0.4704389       NaN
ncp = c(NA, 1:9)
rbeta(n, shape1, shape2, ncp) # [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

Example of bad behaviour for rf:

set.seed(313)
n = 10
df1 = 5
df2 = 20
ncp = c(1:9, NA)
rf(n, df1, df2, ncp) # [1]  0.176153  1.079886  2.086806  2.961054  1.555122 12.506973  1.872629  3.007218  1.297086       NaN
ncp = c(NA, 1:9)
rf(n, df1, df2, ncp) #[1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

The same problem affected stats::rt. It has already been handled at:https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=17306

It appears that stats:rf can be fixed in exactly the same way as stats::rt:

rf_new = function (n, shape1, shape2, ncp)
{
  if (missing(ncp))
    .Call(C_rf, n, df1, df2)
  else (rchisq(n, df1, ncp = ncp)/df1)/(rchisq(n, df2)/df2)
}

rbeta is documented to handle ncp = 0 differently from a missing ncp, so the
default value isn't needed, but I've left it in:

rbeta_new = function (n, shape1, shape2, ncp = 0)
{
  if(missing(ncp)) {
    .Call(C_rbeta, n, shape1, shape2)
  }
  else {
    X <- rchisq(n, 2 * shape1, ncp = ncp)
    X/(X + rchisq(n, 2 * shape2))
  }
}
Comment 1 Martin Maechler 2018-01-11 10:58:28 UTC
Yes, this is a bug, thank you Duncan and Jonas Moss.

--- originally ncp in most of these functions was meant to be of length one (when not missing) and hence  *not* to be recycled.

I will fix this.

-------------------- Related, but not directly on this: -----------

Note that for those of use like me who have the environment variable
	_R_CHECK_LENGTH_1_CONDITION_  set
all this even gives an *error* , already in R 3.4.3:

Note this:

> Sys.unsetenv("_R_CHECK_LENGTH_1_CONDITION_")
> rbeta(2, 5,20, ncp=c(1,NA))
[1] 0.08840147        NaN
Warning messages:
1: In if (is.na(ncp)) { :
  the condition has length > 1 and only the first element will be used
2: In if (ncp == 0) .Call(C_rbeta, n, shape1, shape2) else { :
  the condition has length > 1 and only the first element will be used
3: In rchisq(n, 2 * shape1, ncp = ncp) : NAs produced

> Sys.setenv("_R_CHECK_LENGTH_1_CONDITION_" = TRUE)
> rbeta(2, 5,20, ncp=c(1,NA))
Error in if (is.na(ncp)) { : the condition has length > 1
> 

----->  Can we make _R_CHECK_LENGTH_1_CONDITION = TRUE  the default in R-devel?
Comment 2 Martin Maechler 2018-01-11 16:06:45 UTC
Fixed for R-devel (r 74108) ... to be ported to R-patched in 1-2 days.