[Rd] Error message for infinite probability parameters in rbinom() and rmultinom()

[Christophe Dutang]

Thanks Duncan and Martin for your answer.

> Le 9 avr. 2023 à 00:33, Duncan Murdoch <murdoch.duncan using gmail.com> a écrit :
> 
> On 08/04/2023 5:53 p.m., Martin Maechler wrote:
>>>>>>> Christophe Dutang
>>>>>>>     on Sat, 8 Apr 2023 14:21:53 +0200 writes:
>>     > Dear all,
>>     > Using rmultinom() in a stochastic model, I found this function returns an error message 'NA in probability' for an infinite probability.
>>     > Maybe, a more precise message will be helpful when debugging.
>>     >> rmultinom(1, 3:5, c(1/2, 1/3, Inf))
>>     > Error in rmultinom(1, 3:5, c(1/2, 1/3, Inf)) : NA in probability vector
>>     >> rmultinom(1, 3:5, c(1/2, 1/3, NA))
>>     > Error in rmultinom(1, 3:5, c(1/2, 1/3, NA)) : NA in probability vector
>> Thank you.
>> I agree the first ('Inf') should not do what it currently does,
>> and probably the 2nd one should neither give an error.
>> Note that in rmultinom,  the 'prob' is allowed to be *NOT*
>> scaled to sum(.) = 1.
>> Therefore 'Inf' makes sense as the limit (of a sequence) of (a)
>> very large number(s).
>> I claim that
>>   rmultinom(1, 3, c(1/2, 1/3, Inf))
>> should give the same as
>>   rmultinom(1, 3, c(1/2, 1/3, 1e300))
>> even without a warning,
> 
> That case makes sense, but is it worth the effort?  Certainly
> 
>    rmultinom(1, 3, c(1/2, Inf, Inf))
> 
> can't give a useful answer because we don't know the relative size of the two infinities.  


I think any infinite values in the probability vector or in the size vector should return NaN.


> I imagine the first NA comes from computing prob/sum(prob), which is c(0, 0, NaN).

Line 73 of rmultinom.c check for probability in [0,1] and throw an error by calling C macro ML_WARN_ret_NAN() which return NA or -1.

I also notice that rbinom() throws a warning and not an error (for prob not in [0,1] or missing), e.g. 

> rbinom(4, 2:5, c(1/2, -1, Inf, NA))
[1]  2 NA NA NA

Maybe, only the man page should be updated in the Value field and state that NA is returned when parameter are not correct. By the way, in ?dbinom, it is said 'If size is not an integer, NaN is returned.' which is not true for rbinom() and qbinom(), but only for dbinom() and pbinom().

Thanks in advance for your work

Kind regards, Christophe 

> 
> Duncan Murdoch
> 
> > and OTOH,  an NA in prob may return NA (and signal a warning)
> > instead of an error.
>>     > For rgeom() or rbinom(), we got a warning for infinite probability :
>> Yes, but there, prob must be in [0,1] ... so that's somewhat differnt.
>>     >> rbinom(1, 3, Inf)
>>     > [1] NA
>>     > Warning message:
>>     > In rbinom(1, 3, Inf) : NAs produced
>>     >> rbinom(1, 3, NA)
>>     > [1] NA
>>     > Warning message:
>>     > In rbinom(1, 3, NA) : NAs produced
>>     >> rgeom(1, Inf)
>>     > [1] NA
>>     > Warning message:
>>     > In rgeom(1, Inf) : NAs produced
>>     >> rgeom(1, NA)
>>     > [1] NA
>>     > Warning message:
>>     > In rgeom(1, NA) : NAs produced
>>     > Maybe, it could be better to harmonize the behavior for infinite probability.
>>     > Kind regards, Christophe
>>     >> sessionInfo()
>>     > R version 4.2.3 (2023-03-15)
>>     > Platform: aarch64-apple-darwin20 (64-bit)
>>     > Running under: macOS Ventura 13.2.1
>>     > Matrix products: default
>>     > BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
>>     > LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
>>     > locale:
>>     > [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>>     > attached base packages:
>>     > [1] stats     graphics  grDevices utils     datasets  methods   base
>>     > loaded via a namespace (and not attached):
>>     > [1] compiler_4.2.3 tools_4.2.3
>>     > -------------------------------------------------
>>     > Christophe DUTANG
>>     > LJK, Ensimag, Grenoble INP, UGA, France
>>     > Web: http://dutangc.free.fr
>>     > ______________________________________________
>>     > R-devel using r-project.org mailing list
>>     > https://stat.ethz.ch/mailman/listinfo/r-devel
>> ______________________________________________
>> R-devel using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>