[R] Matrix calculations in R--erroneous?

[Peter Muhlberger]

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 at 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
> 
>