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: ms.unimelb.edu.au/~gabraham