diff --git a/cara/montecarlo.py b/cara/montecarlo.py
index 0a12f297..b7927721 100644
--- a/cara/montecarlo.py
+++ b/cara/montecarlo.py
@@ -86,20 +86,20 @@ def exposure_model_from_vl_breathing_cn():
     fig = plt.figure()
     ax = fig.add_subplot(1, 1, 1)
 
-    n_lines = 30
+    n_lines = 1
     cns = np.linspace(0.01, 0.5, n_lines)
 
     cmap = define_colormap(cns)
 
     for cn in tqdm(cns):
-        er_means = []
+        er_means = np.array([])
         for vl in viral_loads:
             exposure_mc = breathing_exposure_vl(vl)
             exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
             # divide by 2 to have in 30min (half an hour)
             emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
                 cn_B=cn, cn_L=0.2) / 2
-            er_means.append(np.mean(emission_rate))
+            er_means = np.append(er_means, np.mean(emission_rate))
 
         # divide by 2 to have in 30min (half an hour)
         coleman_etal_er_breathing_2 = [x/2 for x in coleman_etal_er_breathing]
@@ -109,15 +109,15 @@ def exposure_model_from_vl_breathing_cn():
             cn, alpha=0.75), linewidth=0.5)
 
     # The dashed line for the chosen Cn,B
-    er_means = []
+    er_means = np.array([])
     for vl in viral_loads:
         exposure_mc = breathing_exposure_vl(vl)
         exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
         # divide by 2 to have in 30min (half an hour)
         emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
             cn_B=0.06, cn_L=0.2) / 2
-        er_means.append(np.mean(emission_rate))
-
+        er_means = np.append(er_means, np.mean(emission_rate))
+    print_er_info(er_means)
     ax.plot(viral_loads, er_means, color=cmap.to_rgba(
         cn, alpha=0.75), linewidth=1, ls='--')
     plt.text(viral_loads[int(len(viral_loads)*0.9)], 10**4.2,
@@ -200,19 +200,19 @@ def exposure_model_from_vl_talking_cn():
     fig = plt.figure()
     ax = fig.add_subplot(1, 1, 1)
 
-    n_lines = 30
+    n_lines = 1
     cns = np.linspace(0.01, 2, n_lines)
     cmap = define_colormap(cns)
 
     for cn in tqdm(cns):
-        er_means = []
+        er_means = np.array([])
         for vl in viral_loads:
             exposure_mc = talking_exposure_vl(vl)
             exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
             # divide by 4 to have in 15min (quarter of an hour)
             emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
                 cn_B=0.1, cn_L=cn) / 4
-            er_means.append(np.mean(emission_rate))
+            er_means = np.append(er_means, np.mean(emission_rate))
 
         # divide by 4 to have in 15min (quarter of an hour)
         coleman_etal_er_talking_2 = [x/4 for x in coleman_etal_er_talking]
@@ -220,16 +220,18 @@ def exposure_model_from_vl_talking_cn():
             cn, alpha=0.75), linewidth=0.5)
 
     # The dashed line for the chosen Cn,L
-    er_means = []
+    er_means = np.array([])
     for vl in viral_loads:
         exposure_mc = talking_exposure_vl(vl)
         exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
         # divide by 4 to have in 15min
         emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
             cn_B=0.06, cn_L=0.2) / 4
-        er_means.append(np.mean(emission_rate))
+        er_means = np.append(er_means, np.mean(emission_rate))
     ax.plot(viral_loads, er_means, color=cmap.to_rgba(
         cn, alpha=0.75), linewidth=1, ls='--')
+    '''Statistical Data'''
+    print_er_info(er_means)
     plt.text(viral_loads[int(len(viral_loads)*0.93)], 10**5.5,
              r"$\mathbf{c_{n,L}=0.2}$", color=cmap.to_rgba(cn), size='small')
 
@@ -370,3 +372,20 @@ def build_breathing_legend(fig):
             width = item.get_window_extent(fig.canvas.get_renderer()).width
             label.set_ha('left')
             label.set_position((-3*width, 0))
+
+
+def print_er_info(er: np.array):
+    """
+    Prints statistical parameters of a given distribution of ER-values
+    :param er: A numpy-array of the ER-values
+    :return: Nothing, parameters are printed
+    """
+    print(f"MEAN of ER = {np.mean(er)}\n"
+          f"SD of ER = {np.std(er)}\n"
+          f"Median of qR = {np.quantile(er, 0.5)}\n")
+ 
+    print(f"Percentiles of ER:")
+    for quantile in (0.01, 0.05, 0.25, 0.50, 0.55, 0.65, 0.75, 0.95, 0.99):
+        print(f"ER_{quantile} = {np.quantile(er, quantile)}")
+
+    return