Re: [R] Matrix calculations in R--erroneous?

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Sat 08 Oct 2005 - 07:55:59 EST

          Rather than adding 1e-15 to all numbers, I suggest you simply make that the floor. (Or use .Machine$double.eps or 2*.Machine$double.eps in place of 1e-15.)

          Another alternative that may or may not apply in your case is to develop an asymptotic expansion for the log(likelihood) for the small numbers. I've had good success with this kind of method. For example, consider the Box-Cox transformation:

          bc(y, b) = (y^b-1)/b

          What do we do with b = 0? We can test for b = 0 and replace those cases by the limit log(y). However, it is numerically more stable to use the following:

          bc(y, b) = ifelse(abs(b*log(y))>.Machine$double.eps, (expm1(b*log(y))/b, log(y)).

          I don't have time to study your example to see if I could see anything like this that could be done, but I think there should be a good chance of finding something like this. Of course, if there are only very few 0's, then it hardly matters. However, if there are quite a few, then you need something like this.

	  hope this helps.
	  spencer graves

Peter Muhlberger wrote:

> Hi Thomas:  Thanks!  Yes, the function
> (yMat2[5,9]-yMat1[5,9]*logitShape()[5,9]) appears to be producing a value of
> -1.102216e-16 rather than 0.  I would have thought it would approach 0 from
> above given that all input values are at or above zero, but evidently not.
> 
> The max function won't do the trick because I need the entire matrix.  I
> could do one cell at a time, but this is part of a ML routine that needs to
> be evaluated hundreds of thousands of times, so I can't afford to slow it
> down that much.
> 
> I guess I can add 1e-15 rather than e-323, but wonder what that might end up
> doing to my estimates.  Guess I'll find out.
> 
> Cheers,  Peter
> 
> On 10/7/05 1:12 PM, "Thomas Lumley" <tlumley@u.washington.edu> wrote:
> 
> 

>>On Fri, 7 Oct 2005, Peter Muhlberger wrote:
>>
>>
>>>Does anyone know how -log(x) can equal 743 but -log(x+0)=Inf? That's what
>>>the following stream of calculations suggest:
>>>
>>>Browse[2]> -log ( 1e-323+yMat2 - yMat1 * logitShape(matrix(parsList$Xs,
>>>nrow = numXs, ncol=numOfCurves), matrix(means, nrow = numXs,
>>>ncol=numOfCurves, byrow=TRUE), matrix(sigmas, nrow = numXs,
>>>ncol=numOfCurves, byrow=TRUE)) )[5,9]
>>>[1] Inf
>>>
>>>Yet:
>>>
>>>Browse[2]> logitShape(matrix(parsList$Xs, nrow = numXs, ncol=numOfCurves),
>>>matrix(means, nrow = numXs, ncol=numOfCurves, byrow=TRUE), matrix(sigmas,
>>>nrow = numXs, ncol=numOfCurves, byrow=TRUE))[5,9]
>>>[1] 1
>>>
>>>So, the logitShape component equals 1.
>>
>>to within 2e-16
>>
>>
>>>Browse[2]> yMat1[5,9]
>>>[1] 1
>>>
>>>So yMat1[5,9]*logitShape()[5,9]=1
>>
>>to within 2e-16
>>
>>
>>>Browse[2]> yMat2[5,9]
>>>[1] 1
>>
>>to within 2e-16
>>
>>
>>>So, yMat2[5,9]-yMat1[5,9]*logitShape()[5,9]=0
>>
>>to within a few parts in 10^16
>>
>>You haven't actually shown us yMat2[5,9]-yMat1[5,9]*logitShape()[5,9],
>>though
>>
>>
>>>Browse[2]> -log ( 1e-323)
>>>[1] 743.7469
>>>
>>>So, -log( 1e-323)=743 while -log( 1e-323+0)=Inf ?
>>>
>>
>>If "0" is really of the order of 1e-16 then this isn't surprising. If the
>>only point of 1e-323 is as a guard value for 0 then use max(1e-323,
>>yMat2[5,9]-yMat1[5,9]*logitShape()[5,9])
>>
>>
>>-thomas
>>
>>
> 
> 
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves@pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Sat Oct 08 08:07:46 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:40:39 EST