qcauchy accuracy (PR#7902)

Testing the code that Morten Welinder suggested for improving
extreme tail behavior of  qcauchy(),
I found what you can read in the subject.
namely that the tan() + floating-point implementation on all
four different versions of Redhat linux, I have access to on
i686 and amd64 architectures,

    > 1/tan(c(-0,0))
gives
    -Inf  Inf

and of course, that can well be considered a feature, since
after all, the tan() function does jump from -Inf to +Inf at 0. 
I was still surprised that this even happens on the R level,
and I wonder if this distinction of "-0" and "0" shouldn't be
mentioned in some place(s) of the R documentation.

For the real problem, the R source (in C), It's simple
to work around the fact that
    qcauchy(0, log=TRUE)
for Morten's code proposal gives -Inf instead of +Inf.

Martin

"MM" == Martin Maechler <maechler@stat.math.ethz.ch>
    on Wed,  1 Jun 2005 08:57:18 +0200 (CEST) writes:

"Morten" == Morten Welinder <mwelinder@gmail.com>
    on Fri, 27 May 2005 20:24:36 +0200 (CEST) writes:

.............

    Morten> Now that pcauchy has been fixed, it is becoming
    Morten> clear that qcauchy suffers from the same problems.

    Morten> 
    Morten> qcauchy(pcauchy(1e100,0,1,FALSE,TRUE),0,1,FALSE,TRUE)

    Morten> should yield 1e100 back, but I get 1.633178e+16.
    Morten> The code below does much better.  Notes:

    Morten> 1. p need not be finite.  -Inf is ok in the log_p
    Morten> case and R_Q_P01_check already checks things.

    MM> yes

    Morten> 2. No need to disallow scale=0 and infinite
    Morten> location.

    MM> yes

    Morten> 3. The code below uses isnan and finite directly.
    Morten> It needs to be adapted to the R way of doing that.

    MM> I've done this, and started testing the new code; a version will
    MM> be put into the next version of R.

    MM> Thank you for the suggestions.

    >>> double
    >>> qcauchy (double p, double location, double scale, int lower_tail, int log_p)
    >>> {
    >>> if (isnan(p) || isnan(location) || isnan(scale))
    >>> return p + location + scale;

    >>> R_Q_P01_check(p);
    >>> if (scale < 0 || !finite(scale)) ML_ERR_return_NAN;

    >>> if (log_p) {
    >>> if (p > -1)
    >>> lower_tail = !lower_tail, p = -expm1 (p);
    >>> else
    >>> p = exp (p);
    >>> }
    >>> if (lower_tail) scale = -scale;
    >>> return location + scale / tan(M_PI * p);
    >>> }

Testing the code that Morten Welinder suggested for improving
extreme tail behavior of  qcauchy(),
I found what you can read in the subject.
namely that the tan() + floating-point implementation on all
four different versions of Redhat linux, I have access to on
i686 and amd64 architectures,

    > 1/tan(c(-0,0))

gives
    -Inf  Inf

and of course, that can well be considered a feature, since
after all, the tan() function does jump from -Inf to +Inf at 0. 
I was still surprised that this even happens on the R level,
and I wonder if this distinction of "-0" and "0" shouldn't be
mentioned in some place(s) of the R documentation.

Indeed I would myself consider this a very useful feature!

However, a query: Clearly from the above (ahich I can reproduce
too), tan() can distinguish between -0 and +0, and return
different results (otherwise 1/tan() would not return different
results).

But how can the user tell the difference between +0 amnd -0?
I've tried the following:

  > sign(c(-0,0))
  [1] 0 0
  > sign(tan(c(-0,0)))
  [1] 0 0
  > sign(1/tan(c(-0,0)))
  [1] -1  1

so sign() is not going to tell us. Is there a function which can?

Short of wrting one's own:

  sign0 <-
    function(x){
      if(abs(x)>0) stop("For this test x must be +0 or -0")
      return(sign(1/tan(x)))
    }

;)

Best wishes,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding@nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 01-Jun-05                                       Time: 10:50:06
------------------------------ XFMail ------------------------------

For the real problem, the R source (in C), It's simple
to work around the fact that
    qcauchy(0, log=TRUE)
for Morten's code proposal gives -Inf instead of +Inf.

Ouch.  Good catch.

Here is what happened: I reduced 1-exp(x) to -expm1(x) which is actual wrong for
x=0 because the results will be differently signed zeros.

I should have reduced 1-exp(x) to 0-expm1(x).

In light of this, any use of R_DT_qIv and R_DT_CIv might have the same problem.
(But, quite frankly, all uses of R_DT_qIv should be eliminated anyway.
 Only qnorm
seems to be using it without killing the right-tail and the log cases.)

Morten

However, a query: Clearly from the above (ahich I can reproduce  
too), tan() can distinguish between -0 and +0, and return different  
results (otherwise 1/tan() would not return different results).

But how can the user tell the difference between +0 amnd -0?

That's indeed a good question - by definition (-0)==(+0) is true,  
-0<0 is false and signum of both -0 and 0 is 0.

I don't see an obvious way of distinguishing them at R level. Besides  
computational ways (like the 1/tan trick) the only (very ugly) way  
coming to my mind is something like:
a==0 && substr(sprintf("%f",a),1,1)=="-"
Note that print doesn't display the sign, only printf does.

At C level it's better - you can use the signbit() function/macro  
there. Any other ideas?

Cheers,
Simon

On Jun 1, 2005, at 5:50 AM, (Ted Harding) wrote:

However, a query: Clearly from the above (ahich I can reproduce
too), tan() can distinguish between -0 and +0, and return different
results (otherwise 1/tan() would not return different results).

But how can the user tell the difference between +0 amnd -0?

That's indeed a good question - by definition (-0)==(+0) is true,
-0<0 is false and signum of both -0 and 0 is 0.

I don't see an obvious way of distinguishing them at R level. Besides
computational ways (like the 1/tan trick) the only (very ugly) way
coming to my mind is something like:
a==0 && substr(sprintf("%f",a),1,1)=="-"
Note that print doesn't display the sign, only printf does.

On my XP machine running R 2.1.0 patched 2005-05-14

sprintf("%f",-0)

[1] "0.000000"

does not print the sign.

however, the tan trick can be done without tan using just division:

R> sign0 <- function(x) if (x != 0) stop("x not zero") else sign(1/x)
R> sign0(0)
[1] 1
R> sign0(-0)
[1] -1