Finished exercise 6.

.vscode/settings.json

@@ -11,7 +11,9 @@
     "ggsave",
     "Khinchine",
     "kleinrock",
+    "linetype",
     "logseq",
+    "maxs",
     "packetsize",
     "pnet",
     "Pollaczek",
@@ -23,6 +25,8 @@
     "typst",
     "xlab",
     "ylab",
-    "ylim"
+    "ylim",
+    "ymax",
+    "ymin"
   ]
 }

code/6.R

@@ -1,3 +1,4 @@
+library(scales)
 library(ggplot2)
 
 calc_avg_delay <- function(ArrivalRate, ServiceRate, StoppingCondition) {
@@ -39,27 +40,57 @@ calc_avg_delay <- function(ArrivalRate, ServiceRate, StoppingCondition) {
   AccumDelay / NumQueueCompleted
 }
 
-generate_graph <- function(rho, output_file) {
-  print(rho)
+generate_graph <- function(ρ_all, output_file) {
+  data <- lapply(ρ_all, function(ρ) {
+    cat(sprintf("ρ = %.2f\n", ρ))
 
-  StoppingConditions <- seq(1000, 2000, by = 500)
-  AvgAvgDelays <- sapply(StoppingConditions, function(StoppingCondition) {
-    sum(sapply(1:50, function(...) {
-      calc_avg_delay(rho, 1, StoppingCondition)
-    }))
+    stopping_conditions <- c(1000, 2000, 5000, 10000, 20000)
+    avg_delay_results <- lapply(stopping_conditions, function(stopping_condition) {
+      cat(sprintf("stopping condition = %.2f\n", stopping_condition))
+      avg_delays <- sapply(1:50, function(...) {
+        calc_avg_delay(ρ, 1, stopping_condition)
+      })
+
+      avg_delays_mean <- mean(avg_delays)
+      avg_delays_var <- var(avg_delays)
+
+      confidence <- 0.95
+      p <- 1 - (1 - confidence) / 2
+      t <- qt(p, length(avg_delays) - 1)
+      avg_delays_ci_min <- avg_delays_mean - t * sqrt(avg_delays_var / length(avg_delays))
+      avg_delays_ci_max <- avg_delays_mean + t * sqrt(avg_delays_var / length(avg_delays))
+
+      list(
+        mean = avg_delays_mean,
+        ci_min = avg_delays_ci_min,
+        ci_max = avg_delays_ci_max
+      )
+    })
+    avg_delay_means <- sapply(avg_delay_results, function(res) res$mean)
+    avg_delay_ci_mins <- sapply(avg_delay_results, function(res) res$ci_min)
+    avg_delay_ci_maxs <- sapply(avg_delay_results, function(res) res$ci_max)
+
+    data.frame(
+      ρ = sprintf("%.2f", ρ),
+      x = stopping_conditions,
+      mean = avg_delay_means,
+      ci_min = avg_delay_ci_mins,
+      ci_max = avg_delay_ci_maxs
+    )
   })
+  data <- do.call("rbind", data)
 
-  data <- data.frame(x = StoppingConditions, y = AvgAvgDelays)
-
-  plot <- ggplot(data) +
-    geom_line(aes(.data$x, .data$y)) +
+  plot <- ggplot(data, aes(.data$x, group = .data$ρ, color = .data$ρ)) +
+    geom_ribbon(aes(ymin = .data$ci_min, ymax = .data$ci_max), alpha = 0.2, fill = "grey50", linetype = 0) +
+    geom_line(aes(y = .data$mean)) +
     xlab("Stopping condition") +
-    ylab("Avg delay")
+    ylab("Avg delay") +
+    scale_x_continuous() +
+    scale_y_log10(labels = scales::comma)
 
   ggsave(plot, file = output_file, device = "svg")
 }
 
 set.seed(0)
 pdf(NULL)
-generate_graph(1, "output/6a.svg")
-generate_graph(2, "output/6b.svg")
+generate_graph(c(0.5, 1, 2), "output/6.svg")

typst/exercises/6.typ

@@ -1,3 +1,22 @@
+#import "/typst/util.typ" as util: indent_par, code_figure
 
+#indent_par[Figure 13 shows the average delay as a function of the stopping condition for system loads $ρ$ equal to 0.5, 1 and 2.]
 
-// TODO: Do a run with rho = 0.5 to compare with.
+#indent_par[We've chosen to include $ρ = 0.5$ to serve as a comparison against the other values]
+
+#figure(
+  image("/output/6.svg", width: 50%),
+	caption: "Average delay as a function of the stopping condition"
+)
+
+#indent_par[The y axis is on a log-scale, because the function grows very rapidly, for certain values we'll discuss shortly. There also exists a ribbon around each value, representing the confidence intervals.]
+
+#indent_par[To calculate the confidence intervals, we perform 50 rounds of the simulator.]
+
+#indent_par[For $ρ = 0.5$, the average delay is constant as the stopping condition grows, indicating this is a stable system.]
+
+#indent_par[For $ρ = 1$ and $ρ = 2$, the average delay increases without bound as the stopping condition grows. This indicates these systems are unstable, with $ρ = 2$ being more unstable than $ρ = 1$.]
+
+#indent_par[As a fun note, another difference between $ρ = 1$ and $ρ = 2$ is that the former has a huge confidence interval, while the latter's is very tight around the mean. This implies that as the system becomes more unstable, the run-to-run stability increases.]
+
+#pagebreak()

typst/main.typ

@@ -62,3 +62,8 @@
 
 === 5. Exercise 5
 #include "exercises/5.typ"
+
+== C. Queueing systems
+
+=== 6. Exercise 6
+#include "exercises/6.typ"

zbuild.yaml

@@ -11,6 +11,7 @@ default:
   - rule: ex5_solve
   - rule: ex5_sim1
   - rule: ex5_sim2
+  - rule: ex6
 
 rules:
   # Typst
@@ -157,3 +158,13 @@ rules:
     exec:
       - - Rscript
         - code/5-sim2.R
+
+  # Exercise 6
+  ex6:
+    out:
+      - output/6.svg
+    deps:
+      - code/6.R
+    exec:
+      - - Rscript
+        - code/6.R