[R] expm() within the Matrix package

[Gad Abraham]

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()

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

(after running your code)

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


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]
[1] 134.5565



-- 
Gad Abraham
Department of Mathematics and Statistics
The University of Melbourne
Parkville 3010, Victoria, Australia
email: g.abraham at ms.unimelb.edu.au
web: http://www.ms.unimelb.edu.au/~gabraham