[R] Matrix multiplication

[peter dalgaard]

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