I would like to find the limiting distribution of my transition matrix after taking n = 30 steps However, my problem was I kept could not get the exact same theoretical value as calculated by hand. My value constantly changes each time I run the code albeit it deviates a little from the theoretical value. Below is the transition matrix that I have configured:
Transition_A
[,1] [,2] [,3]
[1,] 0.29400705 0.7059929 0.0000000
[2,] 0.29400705 0.0000000 0.7059929
[3,] 0.04835626 0.2456508 0.7059929
Now I'm going to run that matrix through a simulation of 1000 trials with n = 30 steps
n<-30 # number of steps to simulate
trials<-1000
#Vectors to store the proportion of simulations in each state
levels_A_G1<-matrix(0, nrow=3, ncol=n+1)
for(j in 1:trials){
X<-1 # initial state of the Markov Chain
levels_A_G1[X, 1] = levels_A_G1[X, 1] + 1 #Record the initial state
for( i in 1:n){
p<-Transition_A[X[i],] # Get the p values
X[i+1]<-sample(x=3, size=1, prob=p) # update the chain
levels_A_G1[X[i+1], i+1] = levels_A_G1[X[i+1], i+1] + 1
}
}
levels_A_G1<-levels_A_G1/trials
levels_A_G1
levels_A_G1[,n+1]
Output:
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[,16] [,17] [,18] [,19] [,20]
[1,] 1 0.259 0.304 0.154 0.183 0.147 0.156 0.118 0.146 0.124 0.138 0.131 0.156 0.133
0.144 0.127 0.153 0.148 0.157 0.166
[2,] 0 0.741 0.172 0.356 0.230 0.278 0.225 0.273 0.245 0.259 0.225 0.261 0.236 0.273
0.240 0.260 0.237 0.250 0.255 0.255
[3,] 0 0.000 0.524 0.490 0.587 0.575 0.619 0.609 0.609 0.617 0.637 0.608 0.608 0.594
0.616 0.613 0.610 0.602 0.588 0.579
[,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31]
[1,] 0.159 0.152 0.142 0.149 0.156 0.161 0.148 0.152 0.151 0.139 0.150
[2,] 0.256 0.257 0.262 0.256 0.262 0.243 0.234 0.257 0.226 0.265 0.242
[3,] 0.585 0.591 0.596 0.595 0.582 0.596 0.618 0.591 0.623 0.596 0.608
#limiting distribution which is the stationary distribution
0.150 0.242 0.608 #Unfortunately this value changes everytime I run my code
Expected output(theoretical value):
0.146, 0.251, 0.603
I understand that my limiting distribution always changes due to the random factor coming from the sample function. If someone could correct my code without using a library package/setting seed. I would gladly appreciate it
For simulating 30 transitions starting from state 1, your code is correct. To get the exact state matrix after 30 steps, take the matrix power:
Transition_A <- matrix(c(0.29400705, 0.29400705, 0.04835626, 0.7059929, 0, 0.2456508, 0, 0.7059929, 0.7059929), 3, 3)
Reduce(`%*%`, replicate(30, Transition_A, simplify = FALSE))
#> [,1] [,2] [,3]
#> [1,] 0.1458786 0.2511173 0.6030028
#> [2,] 0.1458786 0.2511174 0.6030028
#> [3,] 0.1458786 0.2511174 0.6030028
with(eigen(Transition_A), vectors %*% diag(values)^30 %*% solve(vectors))
#> [,1] [,2] [,3]
#> [1,] 0.1458786 0.2511173 0.6030028
#> [2,] 0.1458786 0.2511174 0.6030028
#> [3,] 0.1458786 0.2511174 0.6030028
expm::`%^%`(Transition_A, 30)
#> [,1] [,2] [,3]
#> [1,] 0.1458786 0.2511173 0.6030028
#> [2,] 0.1458786 0.2511174 0.6030028
#> [3,] 0.1458786 0.2511174 0.6030028
matrixcalc::matrix.power(Transition_A, 30)
#> [,1] [,2] [,3]
#> [1,] 0.1458786 0.2511173 0.6030028
#> [2,] 0.1458786 0.2511174 0.6030028
#> [3,] 0.1458786 0.2511174 0.6030028
powerplus::Matpow(Transition_A, 30)
#> [,1] [,2] [,3]
#> [1,] 0.1458786 0.2511173 0.6030028
#> [2,] 0.1458786 0.2511174 0.6030028
#> [3,] 0.1458786 0.2511174 0.6030028
Biodem::mtx.exp(Transition_A, 30)
#> [,1] [,2] [,3]
#> [1,] 0.1458786 0.2511173 0.6030028
#> [2,] 0.1458786 0.2511174 0.6030028
#> [3,] 0.1458786 0.2511174 0.6030028
matpow <- function(x, pow) {
if (pow == 0L) {
1L
} else if(pow %% 2L) {
nextpow <- (pow - 1L)/2L
if (nextpow) {
x %*% Recall(x%*%x, nextpow)
} else x
} else {
Recall(x%*%x, pow/2L)
}
}
matpow(Transition_A, 30)
#> [,1] [,2] [,3]
#> [1,] 0.1458786 0.2511173 0.6030028
#> [2,] 0.1458786 0.2511174 0.6030028
#> [3,] 0.1458786 0.2511174 0.6030028