Skip to content

Matrix calculations in R--erroneous?

6 messages · Spencer Graves, Peter Muhlberger, Thomas Lumley

Original
Peter Muhlberger
# 
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.

Browse[2]> yMat1[5,9]
[1] 1

So yMat1[5,9]*logitShape()[5,9]=1

Browse[2]> yMat2[5,9]
[1] 1

So, yMat2[5,9]-yMat1[5,9]*logitShape()[5,9]=0

Browse[2]> -log (   1e-323)
[1] 743.7469

So, -log( 1e-323)=743 while -log( 1e-323+0)=Inf ?

Any idea of a neat way to get around this?  Even if I put in 1e-50 I still
get Inf.  I deliberately included 1e-323 to insure the function didn't go to
infinity.

Peter
Thomas Lumley
#
On Fri, 7 Oct 2005, Peter Muhlberger wrote:

            

      
      
      
    

        
to within 2e-16

            

      
      
      
    

        
to within 2e-16

            

      
      
      
    

        
to within 2e-16

            

      
      
      
    

        
to within a few parts in 10^16

You haven't actually shown us yMat2[5,9]-yMat1[5,9]*logitShape()[5,9], 
though

            

      
      
      
    

        
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

  
  


  


      

    

    

      
      

        


  

        

      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:

        

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Spencer Graves
    

    

      
      #
    

  


  

      
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:

        

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Peter Muhlberger
    

    

      
      #
    

  


  

      
Hi Spencer:  Thanks!  This gives me a number of other ways of thinking about
this problem.  My one concern is that these approaches would also run into
some difficulties with how long it takes to calculate.  I'm interested not
in a single value but a matrix of over 300k values that has to be recomputed
more than a million times.  If I had to apply ifelse or floor to each of
these, it might take too long.  I'll have to see.

Peter

        
On 10/7/05 5:55 PM, "Spencer Graves" <spencer.graves at pdf.com> wrote:

        

            

      
      
      
    

  
  


  


      

    

    
1 day later

    

      
      

        


  

        

      Thomas Lumley
    

    

      
      #
    

  


  

      
On Fri, 7 Oct 2005, Peter Muhlberger wrote:

            

      
      
      
    

        
pmax, then?

 	-thomas