[R] Log plus one transformation in R

[Duncan Murdoch]

On 12/12/2016 12:26 PM, John Sorkin wrote:
> David,
>   
> I did read the help page. All it says is
> log1p(x) computes log(1+x) accurately also for |x| << 1 (and less
> accurately when x is approximately -1).
> This gives me pause. Does it mean that log(x) does not give accurate
> results? If log1p gives more accurate values than log, why is the log
> function not written to use the more accurate computation performed by
> log1p.

I don't think anyone has directly answered these questions.  The problem 
isn't with log(), it's with the representation of floating point numbers 
in R.  The log of 1 + 10^(-100) is very close to 10^(-100), but since 1 
+ 1e-100 evaluates to 1 (we only keep 15 or 16 digits of precision), 
log(1 + 1e-100) will come out as zero.  On the other hand, log1p(1e-100) 
evaluates correctly to 1e-100.

>   I don't believe I can look directly at the code for log and
> log1p,

R is open source, so you could, but those are likely coming from system 
libraries, so it isn't easy to see how the approximations are being done.

Duncan Murdoch
>   so I need to rely on the kindness of others to explain the
> differences between the computations performed by the functions. I guess
> the test I ran, log1p(0.000001)/log(0.000001+1), did not have enough
> precision to demonstrate a difference between the two functions.
> John
>
>
>
> John David Sorkin M.D., Ph.D.
> Professor of Medicine
> Chief, Biostatistics and Informatics
> University of Maryland School of Medicine Division of Gerontology and
> Geriatric Medicine
> Baltimore VA Medical Center
> 10 North Greene Street
> GRECC (BT/18/GR)
> Baltimore, MD 21201-1524
> (Phone) 410-605-7119
> (Fax) 410-605-7913 (Please call phone number above prior to faxing)
> >>> David Winsemius <dwinsemius at comcast.net> 12/12/16 12:05 PM >>>
>
> > On Dec 12, 2016, at 8:53 AM, John Sorkin <jsorkin at grecc.umaryland.edu>
> wrote:
> >
> > At the risk of being flamed . . .
> > What is the difference between log1p(x) and log(x+1)?
> > The two methods appear to give the same results:
> >> log1p(0.000001)/log(0.000001+1)
> > [1] 1
> > John
>
> Read the help page more carefully.
>