Finished exercise 5.

Co-authored-by: Ricardo-Rodrigues <R1cardu7@users.noreply.github.com>

code/5-sim1.R

@@ -0,0 +1,36 @@
+cur_state <- 1
+trans_matrix <- matrix(
+  c(
+    c(0, 2, 0),
+    c(1, 0, 1),
+    c(1, 1, 0)
+  ),
+  nrow = 3,
+  ncol = 3,
+  byrow = TRUE
+)
+
+
+max_time <- 100000
+cur_time <- 0
+states <- list()
+while (cur_time < max_time) {
+  states[[length(states) + 1]] <- cur_state
+
+  trans_rate <- sum(trans_matrix[cur_state, ])
+  wait_time <- rexp(1, rate = trans_rate)
+  cur_time <- cur_time + wait_time
+
+  probs <- cumsum(trans_matrix[cur_state, ]) / sum(trans_matrix[cur_state, ])
+  prob <- runif(1, 0.0, 1.0)
+  next_state <- which(probs >= prob)
+  if (is.vector(next_state)) {
+    next_state <- next_state[1]
+  }
+  cur_state <- next_state
+}
+
+probs <- sapply(1:3, function(state) length(states[states == state]) / length(states))
+probs <- sapply(probs, function(x) round(x, 4))
+print(probs)
+write.table(probs, "output/5-sim1.csv", sep = "\t", col.names = FALSE, row.names = FALSE)

code/5-sim2.R

@@ -0,0 +1,34 @@
+cur_state <- 1
+trans_matrix <- matrix(
+  c(
+    c(0, 2, 0),
+    c(1, 0, 1),
+    c(1, 1, 0)
+  ),
+  nrow = 3,
+  ncol = 3,
+  byrow = TRUE
+)
+
+
+max_time <- 100000
+cur_time <- 0
+states <- list()
+while (cur_time < max_time) {
+  states[[length(states) + 1]] <- cur_state
+
+  trans_rates <- sapply(1:3, function(state_idx) trans_matrix[cur_state, state_idx])
+
+  # Note: We're suppressing warnings here due to the NAs produced with `trans_rate == 0`
+  wait_times <- suppressWarnings(sapply(trans_rates, function(trans_rate) rexp(1, rate = trans_rate)))
+  next_state <- which.min(wait_times)
+
+  cur_time <- cur_time + wait_times[next_state]
+
+  cur_state <- next_state
+}
+
+probs <- sapply(1:3, function(state) length(states[states == state]) / length(states))
+probs <- sapply(probs, function(x) round(x, 4))
+print(probs)
+write.table(probs, "output/5-sim2.csv", sep = "\t", col.names = FALSE, row.names = FALSE)

code/5-solve.R

@@ -0,0 +1,24 @@
+trans_matrix <- matrix(
+  c(
+    c(0, 2, 0),
+    c(1, 0, 1),
+    c(1, 1, 0)
+  ),
+  nrow = 3,
+  ncol = 3,
+  byrow = TRUE
+)
+
+# Build the coefficient matrix and constants
+coef_matrix <- t(trans_matrix) - diag(rowSums(trans_matrix))
+consts <- c(0, 0, 0)
+
+# Substitute one of the rows by the final `pi0 + pi1 + pi2 = 1` equation.
+subst_idx <- 3
+coef_matrix[subst_idx, ] <- c(1, 1, 1)
+consts[subst_idx] <- 1
+
+solution <- solve(coef_matrix, consts)
+solution <- sapply(solution, function(x) round(x, 4))
+print(solution)
+write.table(solution, "output/5-solve.csv", sep = "\t", col.names = FALSE, row.names = FALSE)

typst/exercises/5.typ

@@ -0,0 +1,101 @@
+#import "/typst/util.typ" as util: indent_par, code_figure
+
+#indent_par[Figure 12 represents the 3-DTMC we'll be using for this exercise.]
+
+#figure(
+  image("/images/5-diagram.png", width: 50%),
+	caption: "3-DTMC"
+)
+
+==== (i) Solving the balance equations
+
+#indent_par[The following equations 6, 7, 8 and 9 are our balance equations:]
+
+$ pi_0 λ_00 + pi_1 λ_10 + pi_2 λ_20 = pi_0 λ_00 + pi_0 λ_01 + pi_0 λ_02 $
+$ pi_0 λ_01 + pi_1 λ_11 + pi_2 λ_21 = pi_1 λ_10 + pi_1 λ_11 + pi_1 λ_12 $
+$ pi_0 λ_02 + pi_1 λ_12 + pi_2 λ_22 = pi_2 λ_20 + pi_2 λ_21 + pi_2 λ_22 $
+$ pi_0 + pi_1 + pi_2 = 1 $
+
+#indent_par[Calling $P$ the matrix $mat(λ_00, λ_01, λ_02; λ_10, λ_11, λ_12; λ_20, λ_21, λ_22)$, the 3 first balance equations can be expressed as the following equation 10:]
+
+$ P^T dot mat(sum_(j) λ_0j, 0, 0; 0, sum_(j) λ_1j, 0; 0, 0, sum_(j) λ_2j;) = 0 $
+
+
+#pagebreak()
+
+#indent_par[We can then model this in R using the following code 2:]
+
+#code_figure(
+  text(size: 0.8em, raw(read("/code/4-solve.R"), lang: "R", block: true)),
+  caption: "Code for solving the balance equations",
+)
+
+#indent_par[Like in exercise 4, we need to substitute one of the equations by the 4th balance equation. Doing so, we reach the limiting state probabilities in table 6:]
+
+#let solution = csv("/output/5-solve.csv", delimiter: "\t")
+#figure(
+	pad(1em, text(size: 1.5em, math.mat(
+		gap: 1em,
+		..solution
+	))),
+	kind: table,
+	caption: "Limiting state probabilities"
+)
+
+#pagebreak()
+
+==== (ii). Simulation using view 1
+
+#indent_par[The following code 4 contains our approach to obtain the limiting state probabilities via simulation implementing the 1st view.]
+
+#code_figure(
+  text(size: 0.8em, raw(read("/code/5-sim1.R"), lang: "R", block: true)),
+  caption: "Code using simulation view 1",
+)
+
+#indent_par[View 1 is similar to how a DTMC works, but before each jump, a state will wait an exponentially distributed amount of time, with rate given by it the transition rate out of it.]
+
+#indent_par[After running the code, we ended up with the results in table 7:]
+
+#let solution = csv("/output/5-sim1.csv", delimiter: "\t")
+#figure(
+	pad(1em, text(size: 1.5em, math.mat(
+		gap: 1em,
+		..solution
+	))),
+	kind: table,
+	caption: "Limiting state probabilities"
+)
+
+#pagebreak()
+
+==== (iii). Simulation using view 2
+
+#indent_par[The following code 5 contains our approach to obtain the limiting state probabilities via simulation implementing the 2nd view.]
+
+#code_figure(
+  text(size: 0.8em, raw(read("/code/5-sim2.R"), lang: "R", block: true)),
+  caption: "Code using simulation view 2",
+)
+
+#indent_par[View 2 instead elects to have multiple wait times, depending on the transition rate to each other state, then selecting the state with the minimum time out of all generated outputs.]
+
+#indent_par[After running the code, we ended up with the results in table 8:]
+
+#let solution = csv("/output/5-sim2.csv", delimiter: "\t")
+#figure(
+	pad(1em, text(size: 1.5em, math.mat(
+		gap: 1em,
+		..solution
+	))),
+	kind: table,
+	caption: "Limiting state probabilities"
+)
+
+==== (iv). Results
+
+#indent_par[Both of the results obtained in tables 7 and 8 are very close to the theoretical values presented in table 6.]
+
+#indent_par[This is, in part, because we used a high limit for the maximum time of 100000. This limit impacts the accuracy of the results greatly, with low maximum times have limiting state probabilities that are very far from the theoretical values.]
+
+#pagebreak()

typst/main.typ

@@ -59,3 +59,6 @@
 
 === 4. Exercise 4
 #include "exercises/4.typ"
+
+=== 5. Exercise 5
+#include "exercises/5.typ"

zbuild.yaml

@@ -8,6 +8,9 @@ default:
   - rule: ex4_solve
   - rule: ex4_matrix
   - rule: ex4_sim
+  - rule: ex5_solve
+  - rule: ex5_sim1
+  - rule: ex5_sim2
 
 rules:
   # Typst
@@ -113,3 +116,33 @@ rules:
     exec:
       - - Rscript
         - code/4-sim.R
+
+  # Exercise 5 solve
+  ex5_solve:
+    out:
+      - output/5-solve.csv
+    deps:
+      - code/5-solve.R
+    exec:
+      - - Rscript
+        - code/5-solve.R
+
+  # Exercise 5 sim1
+  ex5_sim1:
+    out:
+      - output/5-sim1.csv
+    deps:
+      - code/5-sim1.R
+    exec:
+      - - Rscript
+        - code/5-sim1.R
+
+  # Exercise 5 sim2
+  ex5_sim2:
+    out:
+      - output/5-sim2.csv
+    deps:
+      - code/5-sim2.R
+    exec:
+      - - Rscript
+        - code/5-sim2.R