Skip to content

Matrix multiplication

5 messages · annek, Fg Nu, Peter Dalgaard +2 more

#
Hi, 
I have a transition matrix T for which I want to find the steady state matrix for. This could be approximated by taking T^n , for large n. 

T= [ 0.8797   0.0382   0.0527   0.0008
      0.0212    0.8002   0.0041   0.0143
      0.0981    0.0273   0.8802   0.0527
      0.0010    0.1343   0.0630   0.9322]

According to a text book I have T^200 should have reached the steady state L

L =[0.17458813   0.17458813   0.17458813   0.17458813
      0.05731902   0.05731902   0.05731902   0.05731902
      0.35028624   0.35028624   0.35028624   0.35028624
      0.44160126   0.44160126   0.44160126   0.44160126]

I am addressing the problem using a for loop doing matrix multiplication (guess there might be better ways, please suggest) and find a steady state matrix after n=30. But if I run the code with n=100 or more I get "Inf" for all entities in L. Does anyone know why is that?

The code I use look like this
#------------------------------------
rep<-20

T <- Ttest
for(i in 1:rep){
 print(i)
 T<-T%*%Ttest
 Ttest<-T
}
L<-T
print(L)
#----------------------------------
#
Try this

install.packages("expm")
library(expm)

T = matrix(c( 0.8797 , ? 0.0382 ?, 0.0527 ?, 0.0008,
? ? ? ? ? ? ? 0.0212 , ? 0.8002 ?, 0.0041 ?, 0.0143,
? ? ? ? ? ? ? 0.0981 , ? 0.0273 ?, 0.8802 ?, 0.0527,
? ? ? ? ? ? ? 0.0010 , ? 0.1343 ?, 0.0630 ?, 0.9322),?
? ? ? ? ? ?nrow = 4, ncol = 4,byrow = T)
T %^% 200




----- Original Message -----
From: annek <annek at ifm.liu.se>
To: r-help at r-project.org
Cc: 
Sent: Wednesday, December 12, 2012 12:49 PM
Subject: [R] Matrix multiplication

Hi, 
I have a transition matrix T for which I want to find the steady state matrix for. This could be approximated by taking T^n , for large n. 

T= [ 0.8797?  0.0382?  0.0527?  0.0008
? ? ? 0.0212? ? 0.8002?  0.0041?  0.0143
? ? ? 0.0981? ? 0.0273?  0.8802?  0.0527
? ? ? 0.0010? ? 0.1343?  0.0630?  0.9322]

According to a text book I have T^200 should have reached the steady state L

L =[0.17458813?  0.17458813?  0.17458813?  0.17458813
? ? ? 0.05731902?  0.05731902?  0.05731902?  0.05731902
? ? ? 0.35028624?  0.35028624?  0.35028624?  0.35028624
? ? ? 0.44160126?  0.44160126?  0.44160126?  0.44160126]

I am addressing the problem using a for loop doing matrix multiplication (guess there might be better ways, please suggest) and find a steady state matrix after n=30. But if I run the code with n=100 or more I get "Inf" for all entities in L. Does anyone know why is that?

The code I use look like this
#------------------------------------
rep<-20

T <- Ttest
for(i in 1:rep){
print(i)
T<-T%*%Ttest
Ttest<-T
}
L<-T
print(L)
#----------------------------------
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
#
On Dec 12, 2012, at 08:19 , annek wrote:

            
Your code is not reproducible, but I notice that at the end of your loop, T==Ttest, so you are squaring a matrix on every iteration. I.e., you are generating T^2, T^4, T^8, so at i=20 you are computing T^1048576 or so. If you go beyond that, I suppose that even a tiny roundoff will cause an eigenvalue slightly bigger than 1 to blow things up.
#
Suggest you look at the %^% operator from the expm package for raising a matrix to a power.

library(expm)
T %^% 200


Hope this is helpful,

Dan 

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204