###################################################################### # R code: P <- matrix( c( 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.2500, 0.5000, 0.0000, 0.2500, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, 0.0000, 0.0000, 0.0625, 0.2500, 0.1250, 0.2500, 0.2500, 0.0625, 0.0000, 0.0000, 0.0000, 0.2500, 0.5000, 0.2500, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000 ), 6, 6, byrow = T ) apply( P, 1, sum ) # [1] 1 1 1 1 1 1 apply( P, 2, sum ) # [1] 1.3125 0.7500 0.1250 1.7500 0.7500 1.3125 print( P.2 <- P %*% P ) # [,1] [,2] [,3] [,4] [,5] [,6] # [1,] 1.000000 0.0000 0.00000 0.0000 0.0000 0.000000 # [2,] 0.390625 0.3125 0.03125 0.1875 0.0625 0.015625 # [3,] 0.062500 0.2500 0.12500 0.2500 0.2500 0.062500 # [4,] 0.140625 0.1875 0.03125 0.3125 0.1875 0.140625 # [5,] 0.015625 0.0625 0.03125 0.1875 0.3125 0.390625 # [6,] 0.000000 0.0000 0.00000 0.0000 0.0000 1.000000 apply( P.2, 2, sum ) # [1] 1.609375 0.812500 0.218750 0.937500 0.812500 1.609375 print( P.3 <- P.2 %*% P ) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.00000000 0.000000 0.0000000 0.000000 0.000000 0.00000000 [2,] 0.48046875 0.203125 0.0234375 0.171875 0.078125 0.04296875 [3,] 0.14062500 0.187500 0.0312500 0.312500 0.187500 0.14062500 [4,] 0.20703125 0.171875 0.0390625 0.203125 0.171875 0.20703125 [5,] 0.04296875 0.078125 0.0234375 0.171875 0.203125 0.48046875 [6,] 0.00000000 0.000000 0.0000000 0.000000 0.000000 1.00000000 apply( P.3, 2, sum ) # [1] 1.8710938 0.6406250 0.1171875 0.8593750 0.6406250 1.8710938 print( P.4 <- P.3 %*% P ) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.00000000 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 [2,] 0.54199219 0.14453125 0.02148438 0.1367188 0.08203125 0.07324219 [3,] 0.20703125 0.17187500 0.03906250 0.2031250 0.17187500 0.20703125 [4,] 0.26269531 0.13671875 0.02539062 0.1757812 0.13671875 0.26269531 [5,] 0.07324219 0.08203125 0.02148438 0.1367188 0.14453125 0.54199219 [6,] 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000 1.00000000 apply( P.4, 2, sum ) # [1] 2.0849609 0.5351562 0.1074219 0.6523438 0.5351562 2.0849609 print( P.5 <- P.4 %*% P ) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.0000000 0.00000000 0.00000000 0.0000000 0.00000000 0.0000000 [2,] 0.5866699 0.10644531 0.01708984 0.1123047 0.07519531 0.1022949 [3,] 0.2626953 0.13671875 0.02539062 0.1757812 0.13671875 0.2626953 [4,] 0.3078613 0.11230469 0.02197266 0.1376953 0.11230469 0.3078613 [5,] 0.1022949 0.07519531 0.01708984 0.1123047 0.10644531 0.5866699 [6,] 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000 1.0000000 apply( P.5, 2, sum ) # [1] 2.25952148 0.43066406 0.08154297 0.53808594 0.43066406 2.25952148 m.step.transition.matrix <- function( P, m ) { P.m <- P for ( i in 2:m ) { P.m <- P.m %*% P } return( P.m ) } m.step.transition.matrix( P, 5 ) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.0000000 0.00000000 0.00000000 0.0000000 0.00000000 0.0000000 [2,] 0.5866699 0.10644531 0.01708984 0.1123047 0.07519531 0.1022949 [3,] 0.2626953 0.13671875 0.02539062 0.1757812 0.13671875 0.2626953 [4,] 0.3078613 0.11230469 0.02197266 0.1376953 0.11230469 0.3078613 [5,] 0.1022949 0.07519531 0.01708984 0.1123047 0.10644531 0.5866699 [6,] 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000 1.0000000 print( P.10 <- m.step.transition.matrix( P, 10 ) ) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.0000000 0.00000000 0.000000000 0.00000000 0.00000000 0.0000000 [2,] 0.6958740 0.03193378 0.006005764 0.03886700 0.03095722 0.1963623 [3,] 0.4176760 0.04804611 0.009180069 0.05937576 0.04804611 0.4176760 [4,] 0.4333985 0.03886700 0.007421970 0.04804707 0.03886700 0.4333985 [5,] 0.1963623 0.03095722 0.006005764 0.03886700 0.03193378 0.6958740 [6,] 0.0000000 0.00000000 0.000000000 0.00000000 0.00000000 1.0000000 apply( P.10, 2, sum ) # [1] 2.74331069 0.14980412 0.02861357 0.18515682 0.14980412 2.74331069 print( P.100 <- m.step.transition.matrix( P, 100 ) ) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00 [2,] 0.75 1.635836e-10 3.124169e-11 2.022004e-10 1.635836e-10 0.25 [3,] 0.50 2.499335e-10 4.773305e-11 3.089348e-10 2.499335e-10 0.50 [4,] 0.50 2.022004e-10 3.861685e-11 2.499335e-10 2.022004e-10 0.50 [5,] 0.25 1.635836e-10 3.124169e-11 2.022004e-10 1.635836e-10 0.75 [6,] 0.00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.00 apply( P.100, 2, sum ) 3.000000e+00 7.793011e-10 1.488333e-10 9.632691e-10 7.793011e-10 3.000000e+00 ###################################################################### m 1 1.3125 0.7500 0.1250 1.7500 0.7500 1.3125 2 1.609375 0.812500 0.218750 0.937500 0.812500 1.609375 3 1.8710938 0.6406250 0.1171875 0.8593750 0.6406250 1.8710938 4 2.0849609 0.5351562 0.1074219 0.6523438 0.5351562 2.0849609 5 2.25952148 0.43066406 0.08154297 0.53808594 0.43066406 2.25952148 . 10 2.74331069 0.14980412 0.02861357 0.18515682 0.14980412 2.74331069 . 100 3.00000000 0.00000000 0.00000000 0.00000000 0.00000000 3.00000000 ###################################################################### eigen.analysis <- eigen( t( P ) ) $values [1] 1.000000 1.000000 0.809017 0.500000 -0.309017 0.250000 $vectors [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 0.55780321 -3.162278e-01 0.009683312 -0.1360828 [2,] 0 0 -0.32552956 6.324555e-01 -0.265481632 0.5443311 [3,] 0 0 -0.06217061 -7.021667e-17 -0.347519968 -0.2721655 [4,] 0 0 -0.40237667 -2.808667e-16 0.859116607 -0.5443311 [5,] 0 0 -0.32552956 -6.324555e-01 -0.265481632 0.5443311 [6,] 0 1 0.55780321 3.162278e-01 0.009683312 -0.1360828 ###################################################################### # Maple code: with( LinearAlgebra ); P := Matrix( [ [ 1, 0, 0, 0, 0, 0 ], [ 1 / 4, 1 / 2, 0, 1 / 4, 0, 0 ], [ 0, 0, 0, 1, 0, 0 ], [ 1 / 16, 1 / 4, 1 / 8, 1 / 4, 1 / 4, 1 / 16 ], [ 0, 0, 0, 1 / 4, 1 / 2, 1 / 4 ], [ 0, 0, 0, 0, 0, 1 ] ] ); v, e := Eigenvectors( Transpose( P ) ); v; Vector[column](6, [1/4, 1/2, 1, 1, (1/4)*sqrt(5)+1/4, 1/4-(1/4)*sqrt(5)]) e[ 1 .. 6, 3 ]; Vector[column](6, [0, 0, 0, 0, 0, 1]) e[ 1 .. 6, 4 ]; Vector[column](6, [1, 0, 0, 0, 0, 0]) assume( p > 0, p < 1 ); v_equilibrium := Matrix( 1, 6, [ p, 0, 0, 0, 0, ( 1 - p ) ] ); v_equilibrium := Matrix(1, 6, [[p, 0, 0, 0, 0, 1-p]]) Multiply( v_equilibrium, P ); Matrix(1, 6, [[p, 0, 0, 0, 0, 1-p]]) ###################################################################### # R code: rw.sim.1 <- function( k, p, theta.start, n.sim, seed ) { set.seed( seed ) theta <- rep( 0, n.sim + 1 ) theta[ 1 ] <- theta.start for ( i in 1:n.sim ) { theta[ i + 1 ] <- move.1( k, p, theta[ i ] ) if ( ( i %% 10000 ) == 0 ) print( i ) } return( table( theta ) / n.sim ) } move.1 <- function( k, p, theta ) { repeat { increment <- sample( x = c( -1, 0, 1 ), size = 1, prob = p ) theta.next <- theta + increment if ( abs( theta.next ) <= k ) { return( theta.next ) break } } } rw.sim.2 <- function( k, p, theta.start, n.sim, seed ) { set.seed( seed ) theta <- rep( 0, n.sim + 1 ) theta[ 1 ] <- theta.start for ( i in 1:n.sim ) { theta[ i + 1 ] <- move.2( k, p, theta[ i ] ) if ( ( i %% 10000 ) == 0 ) print( i ) } return( table( theta ) / n.sim ) } move.2 <- function( k, p, theta ) { increment <- sample( x = c( -1, 0, 1 ), size = 1, prob = p ) if ( ( theta >= - ( k - 1 ) ) & ( theta <= ( k - 1 ) ) ) { theta.next <- theta + increment } else if ( ( theta == k ) & ( ( theta + increment ) <= k ) ) { theta.next <- theta + increment } else if ( ( theta == k ) & ( ( theta + increment ) > k ) ) { theta.next <- - k } else if ( ( theta == - k ) & ( ( theta + increment ) >= - k ) ) { theta.next <- theta + increment } else if ( ( theta == - k ) & ( ( theta + increment ) < - k ) ) { theta.next <- k } return( theta.next ) } k <- 2 p <- c( 1 / 3, 1 / 3, 1 / 3 ) theta.start <- 0 n.sim <- 100000 seed <- 1 print( distribution.1 <- rw.sim.2( k, p, theta.start, n.sim, seed ) ) [1] 10000 [1] 20000 [1] 30000 [1] 40000 [1] 50000 [1] 60000 [1] 70000 [1] 80000 [1] 90000 [1] 100000 theta -2 -1 0 1 2 0.20057 0.19927 0.19872 0.20087 0.20058 ###################################################################### # R code k <- 2 p <- c( 1 / 3, 1 / 3, 1 / 3 ) theta.start <- 0 n.sim <- 1000000 seed <- 1 print( distribution.1 <- rw.sim.1( k, p, theta.start, n.sim, seed ) ) -2 -1 0 1 2 0.153422 0.231460 0.231413 0.230500 0.153206 ###################################################################### # Maple code: solve( { 3 * q_1 + 2 * q_2 = 1, q_1 = ( 3 / 2 ) * q_2 }, { q_1, q_2 } ); { q_1 = 3 / 13, q_2 = 2 / 13 } evalf( [ 2 / 13, 3 / 13, 3 / 13, 3 / 13, 2 / 13 ] ); [.1538461538, .2307692308, .2307692308, .2307692308, .1538461538] ###################################################################### p <- c( 1 / 3, 1 / 3, 1 / 3 ) k <- 2 P <- matrix( c( p[ 2 ], p[ 3 ], 0, 0, p[ 1 ], p[ 1 ], p[ 2 ], p[ 3 ], 0, 0, 0, p[ 1 ], p[ 2 ], p[ 3 ], 0, 0, 0, p[ 1 ], p[ 2 ], p[ 3 ], p[ 3 ], 0, 0, p[ 1 ], p[ 2 ] ), 5, 5, byrow = T ) P [,1] [,2] [,3] [,4] [,5] [1,] 0.3333333 0.3333333 0.0000000 0.0000000 0.3333333 [2,] 0.3333333 0.3333333 0.3333333 0.0000000 0.0000000 [3,] 0.0000000 0.3333333 0.3333333 0.3333333 0.0000000 [4,] 0.0000000 0.0000000 0.3333333 0.3333333 0.3333333 [5,] 0.3333333 0.0000000 0.0000000 0.3333333 0.3333333 eigen( t( P ) ) $values [1] 1.0000000 0.5393447 0.5393447 -0.2060113 -0.2060113 $vectors [,1] [,2] [,3] [,4] [,5] [1,] -0.4472136 0.6324555 0.000000 0.000000 0.6324555 [2,] -0.4472136 0.1954395 -0.601501 0.371748 -0.5116673 [3,] -0.4472136 -0.5116673 -0.371748 -0.601501 0.1954395 [4,] -0.4472136 -0.5116673 0.371748 0.601501 0.1954395 [5,] -0.4472136 0.1954395 0.601501 -0.371748 -0.5116673 equilibrium.vector <- matrix( rep( -0.4472136, 5 ), 1, 5 ) equilibrium.vector %*% P [,1] [,2] [,3] [,4] [,5] [1,] -0.4472136 -0.4472136 -0.4472136 -0.4472136 -0.4472136 rescaled.equilibrium.vector <- matrix( rep( 0.2, 5 ), 1, 5 ) rescaled.equilibrium.vector %*% P [,1] [,2] [,3] [,4] [,5] [1,] 0.2 0.2 0.2 0.2 0.2 #################################################################