[R] expm() within the Matrix package

[Martin Maechler]

>>>>> "Gad" == Gad Abraham <g.abraham at ms.unimelb.edu.au>
>>>>>     on Fri, 16 Mar 2007 09:48:52 +1100 writes:

    Gad> Laura Hill wrote:
    >> Hi
    >> 
    >> Could anybody give me a bit of advice on some code I'm having trouble with?
    >> 
    >> I've been trying to calculate the loglikelihood of a function iterated over
    >> a data of time values and I seem to be experiencing difficulty when I use
    >> the function expm(). Here's an example of what I am trying to do
    >> 
    >> 
    >> y<-c(5,10)      #vector of 2 survival times
    >> 
    >> p<-Matrix(c(1,0),1,2)   #1x2 matrix
    >> 
    >> Q<-Matrix(c(1,2,3,4),2,2)   #2x2 matrix
    >> 
    >> q<-Matrix(c(1,2),2,1)       #2x1 matrix
    >> 
    >> Loglik<-rep(0,length(y))    #creating empty vector of same length as y
    >> 
    >> for(i in 1:length(y)){
    >> 
    >> Loglik[i]<-log((p %*% (expm(Q*y[i]))) %*% q)  #calculating
    >> 
    >> #  Loglikelihood for each y[i]
    >> }
    >> 
    >> The code works perfectly if I use exp(Q*y[i]) but not for expm()

    Gad> You need to do a type conversion here, because you get the loglik as a 
    Gad> 1x1 Matrix, instead of a scalar:

    Gad> (after running your code)

    >> log(p %*% expm(Q * y[i]) %*% q)
    Gad> 1 x 1 Matrix of class "dgeMatrix"
    Gad> [,1]
    Gad> [1,] 134.5565


    Gad> If you convert to numeric, you can then assign it to Loglik:
    >> Loglik[1] <- as.numeric(log(p %*% expm(Q * y[i]) %*% q))
    >> Loglik[1]
    Gad> [1] 134.5565


Hmm, I don't think that's Laura's problem
(and actually I don't know what her problem is) :

Assignment of a 1 x 1 matrix to a vector is not a problem,
and hence the  as.numeric(.) above  really has no effect :

> ll <- 1:2
> (m <- matrix(pi, 1,1))
         [,1]
[1,] 3.141593
> ll[1] <- m
> ll
[1] 3.141593 2.000000
> 

Martin Maechler, ETH Zurich