[R] Probability(Markov chain transition matrix)

[Bernardo Rangel Tura]

At 17:25 30/04/2004, you wrote:
>So we come up with Transition matrix stating probabilities of each state
>
>FromTo NC   0       1     2     3
>NC      0.79  0.21  0      0    0
>0          0.09  0.73 0.18 0    0
>1          0.09  0.51  0   0.40  0
>2          0.09  0.38  0      0   0.55
>3          0.06  0.32  0      0   0.62
>
>
>Thus if one starts with all the accounts having no credit ð0 =(1,0,0,0,0), after one period the distribution of account is ð1=(0.79, 0.21, 0, 0, 0)  
>after subsequent periods, it becomes
>
>Ð2=(0.64, 0.32, 0.04, 0, 0),
>Ð3= (0.540, 0.387, 0.058, 0.015, 0)
>Ð4=(0.468,0.431, 0.070, 0.023, 0.008)
>Ð5=(0.417, 0.460, 0.077, 0.028, 0.018)
>and ++++
>Ð10=(0.315, 0.512, 0.091, 0.036, 0.046)

Hi 


I create a rotine for a problem like this, I hope this is useful for you
------------
#############################################
##            MARKOV CHAIN                 ## 
##                                         ##
## ini- initial state                      ## 
## trans- transition matrix                ##
## n - number of transitions               ##
## f - ini%*%trans                         ##  
## fase - f consolidate                    ##
##                                         ## 
#############################################



#          initial state

ini<-matrix(c(10,0,0,0,0,0,0,0,0),nrow=3,ncol=3)

#          transition matrix
#
#      [,1] [,2] [,3]
# [1,] 0.85  0.1  0.05
# [2,] 0.00  0.7  0.3
# [3,] 0.00  0.0  1.0
#
# In R the command is

trans<-matrix(c(.85,0,0,.1,.7,0,0.05,0.3,1),ncol=3)

markov<-function(ini,trans,n){
  l<-ncol(ini)
  fase<-matrix(0,nrow=l,ncol=l)
  fases<-array(0,dim=c(nrow(ini),ncol(ini),n))
  for (i in 1:n){
    f<-ini%*%trans
    for (w in 1:l){fase[w,w]<-sum(f[,w])}
  ini<-fase
  fases[,,i]<-fase
# print(fase)
# 
# Fase allow calcule de value for transition:
# If sate 1 value 0,95, state 2 value 0,3 and  state 3
# value 0 QALY. 
# I´m calculate de value of any transition if
#  print(sum(fase%*%c(.95,.3,0)))
#
# 

  }
  return(fases)
}

a<-markov(ini,trans,5)

For your case:

ini<-matrix(c(1,rep(0,24)),ncol=5)
trans<-matrix(c(0.79,0.21,0,0,0,0.09,0.73,0.18,0,0,0.09,0.51,0,0.40,0,0.09,0.38,0,0,0.55,0.06,0.32,0,0,0.62),ncol=5)
solv<-markov(ini,trans,10)
> solv[,,10]
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.3173366 0.0000000 0.0000000 0.0000000 0.0000000
[2,] 0.0000000 0.2876792 0.0000000 0.0000000 0.0000000
[3,] 0.0000000 0.0000000 0.2875772 0.0000000 0.0000000
[4,] 0.0000000 0.0000000 0.0000000 0.2886571 0.0000000
[5,] 0.0000000 0.0000000 0.0000000 0.0000000 0.2810951








Bernardo Rangel Tura, MD, MSc
National Institute of Cardiology Laranjeiras
Rio de Janeiro Brazil