diff --git a/EIVPackage/EIVGeneral/coverage_collect.py b/EIVPackage/EIVGeneral/coverage_collect.py
index 934eaac864beb0425447fa3be2b5442abe3a04a3..642315ebb0935322c83e5b6b8686b44011a38434 100644
--- a/EIVPackage/EIVGeneral/coverage_collect.py
+++ b/EIVPackage/EIVGeneral/coverage_collect.py
@@ -7,7 +7,8 @@ import torch.backends.cudnn
 from EIVGeneral.coverage_metrics import epistemic_coverage
 
 def get_coverages(not_averaged_predictions, y,\
-        q_range=np.linspace(0.1,0.9,num=30)):
+        q_range=np.linspace(0.1,0.9,num=30),
+        noisy_y=True):
     """
     Compute the numerical and theoretical coverage for various coverage
     factors, computed from the `q` in `q_range`, using
@@ -20,6 +21,9 @@ def get_coverages(not_averaged_predictions, y,\
     :param y: A `torch.tensor` of the same shape then the second components
     of `not_averaged_predictions`. If the feature dimension is missing, it is added.
     :param q_range: An iterator through floats between 0 and 1.
+    :param noisy_y: Boolean. If True (the default), `y` is treated as noisy and
+    the total uncertainty is considered. If False, `y` is treated as the
+    unnoisy ground truth.
     :returns: numerical_coverage, theoretical_coverage (numpy arrays)
     """
     numerical_coverage_collection = []
@@ -28,7 +32,7 @@ def get_coverages(not_averaged_predictions, y,\
         numerical_coverage, theoretical_coverage =\
                 epistemic_coverage(\
                 not_averaged_predictions=not_averaged_predictions,
-                y=y, q=q)
+                y=y, q=q, noisy_y=noisy_y)
         numerical_coverage_collection.append(float(numerical_coverage))
         theoretical_coverage_collection.append(float(theoretical_coverage))
     return np.array(numerical_coverage_collection),\
@@ -37,7 +41,8 @@ def get_coverages(not_averaged_predictions, y,\
 
 def get_coverage_distribution(net_iterator, dataloader_iterator,
         device, number_of_draws, q_range=np.linspace(0.1,0.9,num=30),
-        number_of_test_samples=10, stack=True):
+        number_of_test_samples=10, stack=True,
+        noisy_y=True):
     """
     Returns the numerical and theoretical coverage for all nets in
     `net_iterator` with dataloader from `dataloader_iterator`
@@ -53,6 +58,9 @@ def get_coverage_distribution(net_iterator, dataloader_iterator,
     dataloaders in dataloader_iterator to be used.
     :param stack: Boolean. If True (default) the results will be stacked along 
     the last dimension. If False a list will be returned.
+    :param noisy_y: Boolean. If True (the default), `y` is treated as noisy and
+    the total uncertainty is considered. If False, `y` is treated as the
+    unnoisy ground truth.
     :returns: num_coverage_collection, th_coverage_collection
     """
     num_coverage_collection, th_coverage_collection = [],[]
@@ -78,7 +86,8 @@ def get_coverage_distribution(net_iterator, dataloader_iterator,
                 y_collection = torch.concat(y_collection, dim=0)
                 numerical_coverage, theoretical_coverage = get_coverages(
                         not_averaged_predictions=not_av_pred_collection,
-                        y=y_collection, q_range=q_range)
+                        y=y_collection, q_range=q_range,
+                        noisy_y=noisy_y)
                 num_coverage_collection.append(numerical_coverage)
                 th_coverage_collection.append(theoretical_coverage)
     if stack:
diff --git a/Experiments/plot_coverage.py b/Experiments/plot_coverage.py
index 2b625eb6e1c9ea34bf3464af4cc4b5e1bf53055f..f9d3f22e81cc77f0533a9b975d811372bd8ea5c9 100644
--- a/Experiments/plot_coverage.py
+++ b/Experiments/plot_coverage.py
@@ -6,6 +6,7 @@ import importlib
 import os
 import json
 
+import numpy as np
 import torch
 import torch.backends.cudnn
 from torch.utils.data import DataLoader
@@ -16,9 +17,12 @@ from EIVTrainingRoutines import train_and_store
 from EIVGeneral.coverage_collect import get_coverage_distribution
 from EIVGeneral.manipulate_datasets import VerticalCut
 
+# coverages to consider
+q_range = np.linspace(0.1, 0.95)
 
 # load hyperparameters from JSON file
-def get_coverages(data, eiv, number_of_draws, use_ground_truth=False):
+def compute_coverages(data, eiv, number_of_draws,
+        use_ground_truth=False, noisy_y=True):
     """
     Create network and dataloader iterators for `data` (short dataname) and
     feed them into `get_coverage_distribution`.
@@ -43,7 +47,6 @@ def get_coverages(data, eiv, number_of_draws, use_ground_truth=False):
     long_dataname = conf_dict["long_dataname"]
     short_dataname = conf_dict["short_dataname"]
 
-    print(f"Plotting coverage for {long_dataname}")
 
     load_data = importlib.import_module(f'EIVData.{long_dataname}').load_data
 
@@ -75,6 +78,7 @@ def get_coverages(data, eiv, number_of_draws, use_ground_truth=False):
     else:
         train_data, _ = load_data()
 
+    print(f"Computing {'EiV' if eiv else 'non-EiV'} coverage for {long_dataname}")
 
     # train_data only used for finding dimensions 
     input_dim = train_data[0][0].numel()
@@ -161,14 +165,62 @@ def get_coverages(data, eiv, number_of_draws, use_ground_truth=False):
         net_iterator=net_iterator(eiv=eiv),
         dataloader_iterator=dataloader_iterator(),
         device=device,
-        number_of_draws=number_of_draws)
+        number_of_draws=number_of_draws,
+        q_range=q_range,
+        noisy_y = not use_ground_truth)
     return numerical_coverage, theoretical_coverage
 
 
 
-#####
-# numerical_coverage, theoretical_coverage =\
-#     get_coverages(data='quadratic', eiv=True, number_of_draws=100,
-#             use_ground_truth=True)
-# print(numerical_coverage)
-# print(theoretical_coverage)
+####
+
+def plot_data(eiv_coverages, noneiv_coverages, color='r',
+        against_theoretical = False):
+    eiv_numerical_coverage, eiv_theoretical_coverage = eiv_coverages
+    noneiv_numerical_coverage, noneiv_theoretical_coverage = noneiv_coverages
+    if eiv_numerical_coverage is not None and eiv_theoretical_coverage is not None:
+        mean_eiv_numerical_coverage = np.mean(eiv_numerical_coverage, axis=-1) 
+        std_eiv_numerical_coverage = np.std(eiv_numerical_coverage, axis=-1)
+        if against_theoretical:
+            x_values = np.mean(eiv_theoretical_coverage, axis=-1)
+        else:
+            x_values = np.array(q_range)
+        plt.plot(x_values, mean_eiv_numerical_coverage,
+                color=color, linestyle='solid')
+        plt.fill_between(x_values,
+                mean_eiv_numerical_coverage - std_eiv_numerical_coverage,
+                mean_eiv_numerical_coverage + std_eiv_numerical_coverage, 
+                color=color, alpha=0.5)
+    if noneiv_numerical_coverage is not None and noneiv_theoretical_coverage is not None:
+        mean_noneiv_numerical_coverage = np.mean(noneiv_numerical_coverage, axis=-1)
+        std_noneiv_numerical_coverage = np.std(noneiv_numerical_coverage, axis=-1)
+        if against_theoretical:
+            x_values = np.mean(noneiv_theoretical_coverage, axis=-1)
+        else:
+            x_values = np.array(q_range)
+        plt.plot(x_values, mean_noneiv_numerical_coverage,
+                color=color, linestyle='dashed')
+        plt.fill_between(x_values,
+                mean_noneiv_numerical_coverage - std_noneiv_numerical_coverage,
+                mean_noneiv_numerical_coverage + std_noneiv_numerical_coverage,
+                color=color, alpha=0.3)
+
+
+plt.figure(1)
+plt.clf()
+datasets = ['energy', 'kin8nm', 'yacht', 'wine']
+colors = ['red', 'blue', 'green', 'purple']
+use_ground_truth = False
+against_theoretical = True
+for data, color in zip(datasets, colors):
+    eiv_coverages = compute_coverages(data=data, eiv=True,
+            number_of_draws=[100,5],
+            use_ground_truth=use_ground_truth)
+    noneiv_coverages = compute_coverages(data=data, eiv=False,
+            number_of_draws=100,
+            use_ground_truth=use_ground_truth)
+    plot_data(eiv_coverages, noneiv_coverages, color=color,
+            against_theoretical=against_theoretical)
+x_diag = np.linspace(0.0, 1.0)
+plt.plot(x_diag, x_diag, color='k', linestyle='dotted' )
+plt.show()