From cea4edd073a1c1e0efc63d597d1e9d13704e358f Mon Sep 17 00:00:00 2001
From: Sam Hadow <sam.hadow@inbox.lv>
Date: Mon, 17 Mar 2025 19:22:48 +0100
Subject: [PATCH] confidence intervals

---
 src/main.py | 79 +++++++++++++++++++++++++++++++++++++++--------------
 1 file changed, 59 insertions(+), 20 deletions(-)

diff --git a/src/main.py b/src/main.py
index 46d93ee..68767a6 100644
--- a/src/main.py
+++ b/src/main.py
@@ -2,6 +2,8 @@
 import random
 from heapq import heappush, heappop
 import matplotlib.pyplot as plt
+from scipy.stats import t
+import numpy as np
 
 class Event:
     def __init__(self, event_type, request):
@@ -131,38 +133,75 @@ class Simulation:
 
 def simulation_wrapper():
     C_values = [1, 2, 3, 6]
-    simulation_time = 10000
+    simulation_time = 1000
+    num_runs = 10
+    min_runs = 5
+    confidence_level = 0.95
 
-    lambda_vals = [l / 100 for l in range(1, 301)]
+    lambda_vals = [l/100 for l in range(1, 301)]  # λ from 0.01 to 3.00
 
-    plt.figure(figsize=(10, 6))
+    plt.figure(figsize=(12, 8))
 
     for c in C_values:
-        lambda_values = []
-        mean_response_times = []
-        print(f"Processing C={c}...")
-        for lambda_val in lambda_vals:
-            try:
-                sim = Simulation(c, lambda_val)
-                sim.run(simulation_time)
+        lambda_points = []
+        means = []
+        ci_lower = []
+        ci_upper = []
+        print(f"\nProcessing C={c}")
 
-                if not sim.response_times:
-                    print(f"No completed requests for C={c}, λ={lambda_val:.2f}")
+        for lambda_val in lambda_vals:
+            run_results = []
+            loss_rates = []
+
+            # run num_runs simulation for each lambda
+            for _ in range(num_runs):
+                try:
+                    sim = Simulation(c, lambda_val)
+                    sim.run(simulation_time)
+
+                    if len(sim.response_times) > 0:
+                        run_mean = sum(sim.response_times)/len(sim.response_times)
+                        run_results.append(run_mean)
+                        loss_rates.append(sim.loss_rate)
+
+                except ValueError: # lossrate too high
                     continue
 
-                mean_rt = sum(sim.response_times) / len(sim.response_times)
-                lambda_values.append(lambda_val)
-                mean_response_times.append(mean_rt)
-                print(f"C={c}, λ={lambda_val:.2f}, Mean RT={mean_rt:.2f}, Loss Rate={sim.loss_rate:.2%}")
-            except ValueError as e:
-                print(f"Stopped at λ={lambda_val:.2f} for C={c}: {str(e)}")
+            # reject if not enough successful run
+            if len(run_results) >= min_runs:
+                # statistics
+                mean_rt = np.mean(run_results)
+                std_dev = np.std(run_results, ddof=1)
+                n = len(run_results)
+
+                # confidence interval
+                t_value = t.ppf((1 + confidence_level)/2, n-1)
+                ci = t_value * std_dev / np.sqrt(n)
+
+                # loss rate
+                mean_loss = np.mean(loss_rates)
+
+                # store results
+                lambda_points.append(lambda_val)
+                means.append(mean_rt)
+                ci_lower.append(mean_rt - ci)
+                ci_upper.append(mean_rt + ci)
+
+                print(f"C={c}, λ={lambda_val:.2f}, Mean RT={mean_rt:.2f} ± {ci:.2f}, Loss Rate={mean_loss:.2%}")
+            elif len(run_results) > 0:
+                print(f"λ={lambda_val:.2f} skipped - only {len(run_results)} successful run(s)")
+                continue
+            else:
+                print(f"Stopped at λ={lambda_val:.2f} - no successful run")
                 break
 
-        plt.plot(lambda_values, mean_response_times, marker='.', linestyle='-', label=f'C={c}')
+
+        plt.plot(lambda_points, means, label=f'C={c}')
+        plt.fill_between(lambda_points, ci_lower, ci_upper, alpha=0.2)
 
     plt.xlabel('Arrival Rate (λ)')
     plt.ylabel('Mean Response Time')
-    plt.title('Mean Response Time depending on Arrival Rate for different configurations of C (number of clusters)')
+    plt.title(f'Mean Response Time vs Arrival Rate ({num_runs} runs, 95% CI)')
     plt.legend()
     plt.grid(True)
     plt.show()