From 5a601c9155a872e230cdc07bc9ee113a6e9029e8 Mon Sep 17 00:00:00 2001
From: Sam Hadow <sam.hadow@inbox.lv>
Date: Fri, 9 May 2025 10:06:07 +0200
Subject: [PATCH] range loss simulation + auto select normal distribution or
 student

---
 src/simulation.py | 30 ++++++++++++++++++++++--------
 1 file changed, 22 insertions(+), 8 deletions(-)

diff --git a/src/simulation.py b/src/simulation.py
index 54e2ecc..e2c33d4 100644
--- a/src/simulation.py
+++ b/src/simulation.py
@@ -3,7 +3,7 @@ import os
 import time
 from heapq import heappush, heappop
 import matplotlib.pyplot as plt
-from scipy.stats import t
+from scipy.stats import t, norm
 import numpy as np
 from multiprocessing import Pool
 
@@ -182,8 +182,13 @@ def simulation_wrapper1(simulation_time, num_runs, min_runs, confidence_level, l
                     std_dev = np.std(run_results, ddof=1)
 
                     # confidence interval
-                    t_value = t.ppf((1 + confidence_level)/2, n-1)
-                    ci = t_value * std_dev / np.sqrt(n)
+                    if n>= 30:
+                        z_value = norm.ppf((1 + confidence_level)/2)
+                        ci = z_value * std_dev / np.sqrt(n)
+                    else:
+                        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)
@@ -245,7 +250,11 @@ def simulation_wrapper2(simulation_time, num_runs, min_runs, confidence_level, l
                 f"{first_lambda if first_lambda is not None else 'never'}; "
                 f" last λ with data: {last_lambda}")
 
-            lambda_vals2.extend([last_lambda + i/1000 for i in range(-100,51)])
+            if first_lambda is not None:
+                range_vals = int((last_lambda - first_lambda + 0.5)*1000)
+                lambda_vals2.extend(first_lambda + i/1000 for i in range(0, range_vals) )
+            else:
+                lambda_vals2.extend([last_lambda + i/1000 for i in range(-100,51)])
 
 
 
@@ -271,8 +280,12 @@ def simulation_wrapper2(simulation_time, num_runs, min_runs, confidence_level, l
                         std_dev = np.std(loss_rates, ddof=1)
 
                         # confidence interval
-                        t_value = t.ppf((1 + confidence_level)/2, n-1)
-                        ci = t_value * std_dev / np.sqrt(n)
+                        if n>= 30:
+                            z_value = norm.ppf((1 + confidence_level)/2)
+                            ci = z_value * std_dev / np.sqrt(n)
+                        else:
+                            t_value = t.ppf((1 + confidence_level)/2, n-1)
+                            ci = t_value * std_dev / np.sqrt(n)
 
                         # store results
                         lambda_points.append(lambda_val)
@@ -297,7 +310,7 @@ def simulation_wrapper2(simulation_time, num_runs, min_runs, confidence_level, l
         sorted_C = sorted(mean_at_lambda1.items(), key=lambda item: item[1][0])
         (best_C, (mean1, lower1, upper1)), (_, (_, lower2, _)) = sorted_C[:2]
         if upper1 < lower2:
-            print(f"Optimal C at λ=1 is {best_C} with Mean RT = {mean1:.2f} (non-overlapping CIs)")
+            print(f"\nOptimal C at λ=1 is {best_C} with Mean RT = {mean1:.2f} (non-overlapping CIs)")
         else:
             print("Confidence intervals overlap between the two best C.")
     else:
@@ -305,10 +318,11 @@ def simulation_wrapper2(simulation_time, num_runs, min_runs, confidence_level, l
 
     #plot curves
     if len(lambda_vals) > 1:
-        ax_loss.set_xlim(left=2)
+        ax_loss.set_xlim(left=1)
         ax_loss.set_xlabel('Arrival Rate (λ)')
         ax_loss.set_ylabel('Mean Loss Rate')
         ax_loss.set_title(f'Mean Loss Rate vs Arrival Rate ({num_runs} runs, {int(confidence_level*100)}% CI)')
+        ax_loss.axhline(y=0.05, color='red', linestyle='--', linewidth=1, label='5% Loss Threshold')
         ax_loss.legend()
         ax_loss.grid(True)