[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



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

    Gad> ______________________________________________
    Gad> R-help at stat.math.ethz.ch mailing list
    Gad> https://stat.ethz.ch/mailman/listinfo/r-help
    Gad> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
    Gad> and provide commented, minimal, self-contained, reproducible code.