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:> 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) > #----------------------------------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. -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
> -----Original Message----- > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r- > project.org] On Behalf Of annek > Sent: Tuesday, December 11, 2012 11:19 PM > To: r-help at r-project.org > 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) andSuggest 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
One solution that does not require matrix multiplication: Remember that the steady state vector is in the nullspace of I - T. Therefore: require(MASS) n1 <- Null(t(diag(nrow(T)) - T)) n1 / sum(n1) On Wed, Dec 12, 2012 at 2:19 AM, annek <annek@ifm.liu.se> wrote:> 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@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. > > >[[alternative HTML version deleted]]