Markov Chain Chap 9 - 21

> trmatrix <- function(n, p, r) {
+     p_l <- r * (1 - p)
+     p_r <- p * (1 - r)
+     p_0 <- r * p + (1 - p) * (1 - r)
+     Q <- matrix(data = 0, nrow = n, ncol = n)
+     diag(Q) <- rep(p_0, n)
+     for (i in 1:(n - 1)) Q[i, (i + 1)] <- p_r
+     for (i in 2:n) Q[i, (i - 1)] <- p_l
+     Q[1, 1] <- 1 - p
+     Q[1, 2] <- p
+     Q[n, n] <- 1 - r
+     Q[n, (n - 1)] <- r
+     return(Q)
+ }

> simula <- function(trans, N) {
+     transita <- function(char, trans) {
+         sample(colnames(trans), 1, prob = trans[char, ])
+     }
+     sim <- character(N)
+     sim[1] <- sample(colnames(trans), 1)
+     for (i in 2:N) {
+         sim[i] <- transita(sim[i - 1], trans)
+     }
+     sim
+ }

> n <- 20
> p <- 0.6
> r <- 0.3
> Q1 <- trmatrix(n, p, r)
> rownames(Q1) <- 1:20
> colnames(Q1) <- 1:20
> Q1Simulate <- as.numeric(simula(Q1, 1000))
> n <- 20
> p <- 0.3
> r <- 0.6
> Q2 <- trmatrix(n, p, r)
> rownames(Q2) <- 1:20
> colnames(Q2) <- 1:20
> Q2Simulate <- as.numeric(simula(Q2, 1000))
> n <- 20
> p <- 0.5
> r <- 0.5
> Q3 <- trmatrix(n, p, r)
> rownames(Q3) <- 1:20
> colnames(Q3) <- 1:20
> Q3Simulate <- as.numeric(simula(Q3, 1000))
> n <- 20
> p <- 0.3
> r <- 0.3
> Q4 <- trmatrix(n, p, r)
> rownames(Q4) <- 1:20
> colnames(Q4) <- 1:20
> Q4Simulate <- as.numeric(simula(Q4, 1000))

> par(mfrow = c(2, 2))
> plot(Q1Simulate, type = "l", col = "blue", main = ("traffic intensity 0.6/0.3"))
> plot(Q2Simulate, type = "l", col = "blue", main = ("traffic intensity 0.3/0.6"))
> plot(Q3Simulate, type = "l", col = "blue", main = ("traffic intensity 0.5/0.5"))
> plot(Q4Simulate, type = "l", col = "blue", main = ("traffic intensity 0.3/0.3"))

> par(mfrow = c(2, 2))
> plot(table(Q1Simulate)/10)
> plot(table(Q2Simulate)/10)
> plot(table(Q3Simulate)/10)
> plot(table(Q4Simulate)/10)

> mean(Q1Simulate)
[1] 19.121
> mean(Q2Simulate)
[1] 2.079
> mean(Q3Simulate)
[1] 15.547
> mean(Q4Simulate)
[1] 12.575