diff --git a/Experiments/plot_coverage.py b/Experiments/plot_coverage.py
index f4a93a4c8776b4101bc6277e885d18e3fdf786b3..298603ce4b400606fa21e16c7997410183b7973f 100644
--- a/Experiments/plot_coverage.py
+++ b/Experiments/plot_coverage.py
@@ -12,6 +12,7 @@ import numpy as np
 import torch
 import torch.backends.cudnn
 from torch.utils.data import DataLoader
+from matplotlib.pyplot import cm
 import matplotlib.pyplot as plt
 
 from EIVArchitectures import Networks
@@ -24,7 +25,7 @@ q_range = np.linspace(0.1, 0.95)
 
 # load hyperparameters from JSON file
 def compute_coverages(data, eiv, number_of_draws,
-        use_ground_truth=False, noisy_y=True):
+        use_ground_truth=False):
     """
     Create network and dataloader iterators for `data` (short dataname) and
     feed them into `get_coverage_distribution`.
@@ -191,32 +192,44 @@ def coverage_diagonal_plot(eiv_coverages, noneiv_coverages, color,
     """
     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', label=f'eiv_{label}')
-        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', label=f'noneiv_{label}')
-        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)
+    assert (len(eiv_numerical_coverage.shape)) == 2
+    assert (len(noneiv_numerical_coverage.shape)) == 2
+    # EiV
+    # take mean/std over seed dimension
+    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:
+        # show theoretical coverage on x-axis
+        x_values = np.mean(eiv_theoretical_coverage, axis=-1)
+    else:
+        # show q-range on x-axis
+        x_values = np.array(q_range)
+    # plot mean
+    plt.plot(x_values, mean_eiv_numerical_coverage,
+            color=color, linestyle='solid', label=f'eiv_{label}')
+    # plot std
+    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)
+    # non-EiV
+    # take mean/std over seed dimension
+    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:
+        # show theoretical coverage on x-axis
+        x_values = np.mean(noneiv_theoretical_coverage, axis=-1)
+    else:
+        # show q-range on x-axis
+        x_values = np.array(q_range)
+    # plot mean
+    plt.plot(x_values, mean_noneiv_numerical_coverage,
+            color=color, linestyle='dashed', label=f'noneiv_{label}')
+    # plot std
+    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)
 
 def coverage_residual_plot(eiv_coverages, noneiv_coverages, color,
         absolute_values=True, against_theoretical=False, label=''):
@@ -275,7 +288,7 @@ def coverage_residual_plot(eiv_coverages, noneiv_coverages, color,
 # create figures, together with title and axis labels
 plt.figure(1)
 plt.clf()
-plt.title('Coverage for datasets with ground truth')
+plt.title('Coverage for datasets with ground truth (vs. q)')
 plt.xlabel('q')
 plt.ylabel('coverage')
 plt.figure(2)
@@ -285,23 +298,24 @@ plt.xlabel('th. coverage')
 plt.ylabel('coverage')
 plt.figure(3)
 plt.clf()
-plt.title('Coverage of noisy labels (vs. q)')
+plt.title('q Deviation of coverage of noisy labels')
 plt.xlabel('q')
 plt.ylabel('coverage')
 plt.figure(4)
 plt.clf()
-plt.title('Theoretical deviation Coverage for datasets with ground truth')
+plt.title('q Deviation of coverage for datasets with ground truth')
 plt.xlabel('q')
 plt.ylabel('deviation cov from q')
 plt.figure(5)
 plt.clf()
-plt.title('Deviation (theory) of noisy Coverage')
+plt.title('Theory deviation  of noisy coverage')
 plt.xlabel('q')
 plt.ylabel('deviation cov from theor. cov')
 
 # datasets to plot and their coloring
-datasets = ['linear', 'quadratic','yacht','wine','power']
-colors = ['red', 'blue','green','purple','orange','cyan']
+datasets = ['linear', 'quadratic','yacht','wine','power',
+        'protein','concrete','california','energy','kin8nm','msd','naval']
+colors = cm.rainbow(np.linspace(0,1,len(datasets)))
 
 # loop through data
 for use_ground_truth in [False, True]:
@@ -354,5 +368,17 @@ for fig_nr in range(1,6):
     plt.figure(fig_nr)
     plt.legend()
 
-# show plots
+plt.figure(1)
+plt.savefig('results/figures/summary_coverage_ground_truth.pdf')
+plt.figure(2)
+plt.savefig('results/figures/summary_coverage_noisy.pdf')
+plt.figure(3)
+plt.savefig('results/figures/summary_q_deviation_coverage_noisy.pdf')
+plt.figure(4)
+plt.savefig('results/figures/summary_q_deviation_coverage_ground_truth.pdf')
+plt.figure(5)
+plt.savefig('results/figures/summary_theory_deviation_coverage_ground_truth.pdf')
+
+# show plot
 plt.show()
+
diff --git a/README.md b/README.md
index 7547a9ddf61f21fe3abacba7584910cedd0dfe34..11bb840651f4a26d5b4fc09d78555994dc92f5b5 100644
--- a/README.md
+++ b/README.md
@@ -36,7 +36,7 @@ pip install EIVPackage/
 Installing this package will make 3 modules available to the python environment: `EIVArchitectures` (for building EiV Models), `EIVTrainingRoutines` (containing a general training framework), `EIVGeneral` (containing a single module needed for repeated sampling).
 
 ## Missing
-+ Tell to create folders, like `Experiments/results`
++ Tell to create folders, like `Experiments/results`,  `Experiments/results/figures/`
 
 ## Contributing