[Rd] NA warnings for r<distr>() {aka "patch for random.c"}

Martin Maechler maechler at stat.math.ethz.ch
Wed Mar 12 08:39:18 CET 2008


>>>>> "BAT" == Berwin A Turlach <statba at nus.edu.sg>
>>>>>     on Wed, 12 Mar 2008 01:23:10 +0800 writes:

    BAT> G'day Martin, On Tue, 11 Mar 2008 18:07:35 +0100 Martin
    BAT> Maechler <maechler at stat.math.ethz.ch> wrote:

    >> >>>>> "BAT" == Berwin A Turlach <statba at nus.edu.sg> >>>>>
    >> on Tue, 11 Mar 2008 13:19:46 +0800 writes:

    BAT> [...]  The first two lines give identical results, as
    BAT> one could reasonably expect.
    >> 
    >> Yes, but I don't think a user should *rely* on the way R
    >> generates these random number; though I agree for the
    >> very specific case of rlnorm(..).
    >> 
    BAT> But the other two do not and I would argue that a user
    BAT> could reasonably expect that the commands in these two
    BAT> lines should lead to identical results.
    >> 
    >> They now do.

    BAT> Well, actually I forgot to mention one could also put
    BAT> another argument forward.  As the log-mean parameter of
    BAT> the lognormal distribution goes to -Inf, the
    BAT> distribution degenerates to something that has mean 0
    BAT> and variance 0, i.e. could be taken as the constant
    BAT> zero and, hence, one might expect that rlnorm(1, -Inf)
    BAT> returns 0.

yes; and that has been a consequence of the changes I did yesterday.

    BAT> But as the log-mean parameter goes to Inf, the
    BAT> distribution degenerates to something with infinite
    BAT> mean and infinite variance.  Thus, perhaps it is more
    BAT> sensible for rlnorm(1, Inf) to return NaN instead of
    BAT> Inf.....

well.  We have pretty strictly tried to follow the principle
that 'Inf' or '-Inf'  should be used when there are clear limit

If     lim_{x \to x_0} f(x) =  +/- Inf =: M

we would set
		f(x_0) = M

for x_0 in   |R + {-Inf, Inf}

If we apply this principle here, clearly,

rlnorm(., Inf) :=  lim_{lmu -> Inf}   rlnorm(., lmu)   =
	        =  lim_{lmu -> Inf} exp(rnorm(., lmu)) =
	        =  exp(lim_{lmu -> Inf} rnorm(., lmu)) =
	        =  exp(Inf) = Inf

{ or you directly numerically see that already
  rlnorm(50, 1000)   gives all Inf }

Martin

    >> I don't think your change to .../R/distn.R was good,

    BAT> I didn't like it either, but it was the simplest way I
    BAT> could think of that would allow the C rexp() routine to
    BAT> realise that a scale parameter of 0 actually came from
    BAT> a rate parameter of -Inf in the R code.

    >> but the others I have more or less committed together
    >> with a few more similar ones.

    BAT> Thanks.

    BAT> BTW, I was surprised to realise that the *exp()
    BAT> functions in the underlying C code use the opposite
    BAT> parameterisation from the corresponding functions at R
    BAT> level.  Perhaps it would be worthwhile to point this
    BAT> out in section 6.7.1 of the Writing R extension manual?
    BAT> In particular since the manual states:
    >> 
    BAT> Note that these argument sequences are (apart from the
    BAT> names and that @code{rnorm} has no @var{n}) exactly the
    BAT> same as the corresponding @R{} functions of the same
    BAT> name, so the documentation of the @R{} functions can be
    BAT> used.
    >> 
    BAT> Well, as I noticed the hard way, for *exp() the
    BAT> documentation of the corresponding R functions cannot
    BAT> be used. ;-)
    >> 
    >> We often also gratefully accept patches for the
    >> documentation

    BAT> I know, and I am always amazed that despite this policy
    BAT> (or perhaps because of it?) the documentation of R is
    BAT> not patchy.... ;-)

    BAT> Cheers,
	
    BAT> 	Berwin



More information about the R-devel mailing list