diff --git a/EIVPackage/EIVGeneral/coverage_collect.py b/EIVPackage/EIVGeneral/coverage_collect.py
new file mode 100644
index 0000000000000000000000000000000000000000..642315ebb0935322c83e5b6b8686b44011a38434
--- /dev/null
+++ b/EIVPackage/EIVGeneral/coverage_collect.py
@@ -0,0 +1,96 @@
+"""
+Collect coverages for various coverage factors, networks and dataloaders.
+"""
+import numpy as np
+import torch
+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),
+        noisy_y=True):
+    """
+    Compute the numerical and theoretical coverage for various coverage
+    factors, computed from the `q` in `q_range`, using
+    `EIVGeneral.epistemic_coverage` 
+    :param not_averaged_predictions: A tuple of tensors as in the output of
+    `FNNEIV.predict` with `take_average_of_prediction` set to `False`, i.e.:
+    the predictions of the neural net not averaged over the first dimension 
+    (the repetition dimension in `FNNEIV.predict`) and
+    the aleatoric uncertainty with a batch dimension and a feature dimension. 
+    :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 = []
+    theoretical_coverage_collection = []
+    for q in q_range:
+        numerical_coverage, theoretical_coverage =\
+                epistemic_coverage(\
+                not_averaged_predictions=not_averaged_predictions,
+                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),\
+            np.array(theoretical_coverage_collection)
+    
+
+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,
+        noisy_y=True):
+    """
+    Returns the numerical and theoretical coverage for all nets in
+    `net_iterator` with dataloader from `dataloader_iterator`
+    :param net_iterator: Iterator yielding `torch.nn.Module`
+    :param dataloader_iterator: Iterator yielding `torch.utils.data.DataLoader`.
+    Should return the same number of elements than `net_iterator`.
+    :param device: The `torch.device` to be used.
+    :param number_of_draws: Integer, to be used as keyword for all networks
+    returned by `net_iterator`. Defaults to 10.
+    :param q_range: Range of floats between 0 and 1, as in 
+    `EIVGeneral.coverage_plot.get_coverages`.
+    :param number_of_test_samples: An integer, number of samples from the
+    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 = [],[]
+    for net, dataloader in\
+     zip(net_iterator, dataloader_iterator):
+                not_av_pred_collection_out, not_av_pred_collection_sigma,\
+                    y_collection = [], [], []
+                for i, (x,y) in enumerate(dataloader):
+                    if i>= number_of_test_samples:
+                        break
+                    x, y = x.to(device), y.to(device)
+                    not_averaged_predictions = net.predict(x,
+                            take_average_of_prediction=False,
+                            number_of_draws=number_of_draws)
+                    not_av_pred_collection_out.append(
+                            not_averaged_predictions[0])
+                    not_av_pred_collection_sigma.append(
+                            not_averaged_predictions[1])
+                    y_collection.append(y)
+                not_av_pred_collection = [
+                        torch.concat(not_av_pred_collection_out, dim=0), 
+                        torch.concat(not_av_pred_collection_sigma, dim=0)]
+                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,
+                        noisy_y=noisy_y)
+                num_coverage_collection.append(numerical_coverage)
+                th_coverage_collection.append(theoretical_coverage)
+    if stack:
+        num_coverage_collection = np.stack(num_coverage_collection, axis=-1)
+        th_coverage_collection = np.stack(th_coverage_collection, axis=-1)
+    return num_coverage_collection, th_coverage_collection
diff --git a/EIVPackage/EIVGeneral/manipulate_datasets.py b/EIVPackage/EIVGeneral/manipulate_datasets.py
new file mode 100644
index 0000000000000000000000000000000000000000..f04280e3a84f374cef1705a3e53a6607a6d0e394
--- /dev/null
+++ b/EIVPackage/EIVGeneral/manipulate_datasets.py
@@ -0,0 +1,21 @@
+from operator import itemgetter
+from torch.utils.data import Dataset
+
+
+class VerticalCut(Dataset):
+    """
+    Only include certain elements of the dataset `dataset`, namely those with
+    index in `components_to_pick`.
+    :param dataset: A torch.utils.data.Dataset
+    :param components_to_pick: List, defaults to [0,].
+    """
+    def __init__(self, dataset, components_to_pick=[0,]):
+        super(VerticalCut, self).__init__()
+        self.dataset = dataset
+        self.components_to_pick = components_to_pick
+
+    def __getitem__(self, i):
+        return itemgetter(*self.components_to_pick)(self.dataset.__getitem__(i))
+
+    def __len__(self):
+        return len(self.dataset)
diff --git a/Experiments/create_tabular.py b/Experiments/create_tabular.py
index 728b6e5228a3e1849210ff5a4552d4a35873de65..8b8373de4c5735c380cde7b9079b39e224a0af53 100644
--- a/Experiments/create_tabular.py
+++ b/Experiments/create_tabular.py
@@ -1,9 +1,8 @@
 import os
 import glob
-import argparse
 import json
 
-metrics_to_display = ['rmse','logdens','bias','true_coverage_numerical','avg_bias']
+metrics_to_display = ['rmse','coverage_theory','coverage_numerical','true_coverage_numerical','avg_bias']
 show_incomplete = True
 
 list_of_result_files = glob.glob(os.path.join('results','*.json'))
diff --git a/Experiments/evaluate_metrics.py b/Experiments/evaluate_metrics.py
index 22d4e66495a8e87cbc6230c825d728629cb53ae2..bf052e872d2e42d830e827526e52649da4d67be4 100644
--- a/Experiments/evaluate_metrics.py
+++ b/Experiments/evaluate_metrics.py
@@ -20,7 +20,7 @@ from EIVData.repeated_sampling import repeated_sampling
 
 # read in data via --data option
 parser = argparse.ArgumentParser()
-parser.add_argument("--data", help="Loads data", default='replin')
+parser.add_argument("--data", help="Loads data", default='linear')
 parser.add_argument("--no-autoindent", help="",
         action="store_true") # to avoid conflics in IPython
 args = parser.parse_args()
@@ -285,7 +285,8 @@ def collect_full_seed_range_metrics(load_data,
                 hidden_layers = noneiv_conf_dict["hidden_layers"]
                 saved_file = os.path.join('saved_networks',
                             f'noneiv_{short_dataname}'\
-                                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                                    f'_init_std_y_{init_std_y:.3f}'\
+                                    f'_ureg_{unscaled_reg:.1f}'\
                                     f'_p_{p:.2f}_seed_{seed}.pkl')
                 net = Networks.FNNBer(p=p, init_std_y=init_std_y,
                         h=[input_dim, *hidden_layers, output_dim]).to(device)
@@ -455,13 +456,7 @@ seed_list = range(noneiv_conf_dict["seed_range"][0],
 
 
 
-
-
-
-
-
-
-max_batch_number = 2
+number_of_test_samples = 2
 for seed in tqdm(seed_list):
     try:
         train_data, test_data, true_train_data, true_test_data \
@@ -478,7 +473,7 @@ for seed in tqdm(seed_list):
                 batch_size=int(np.min((len(true_test_data), 800))), shuffle=True)
     for i in tqdm(range(num_test_epochs)):
         for j, x_y_pairs in enumerate(test_dataloader):
-            if j > max_batch_number:
+            if j > number_of_test_samples:
                 break
             # fill in ground truth with None, if not existent
             if true_test_data is None:
diff --git a/Experiments/plot_coverage.py b/Experiments/plot_coverage.py
new file mode 100644
index 0000000000000000000000000000000000000000..b7a511645c34f2dd62f693e037db36ec05111696
--- /dev/null
+++ b/Experiments/plot_coverage.py
@@ -0,0 +1,442 @@
+"""
+Compute the coverage for various coverage factors and compare them with 
+with the corresponding q or theoretical coverage. Results will be stored
+in various plots in the results/figures folder.
+"""
+import importlib
+import os
+import json
+
+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
+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 compute_coverages(data, eiv, number_of_draws,
+        use_ground_truth=False):
+    """
+    Create network and dataloader iterators for `data` (short dataname) and
+    feed them into `get_coverage_distribution`.
+    :data: String, short dataname
+    :eiv: Boolean. If True an EiV model is used, else an non-EiV model.
+    :number_of_draws: Number of draws to use for prediction. Take an int for
+    non-EiV models and a two-element list for EiV models.
+    :use_ground_truth: If True, unnoisy `y` are considered when computing the
+    coverage. If there is no ground truth None,None is returned
+    :returns: numerical_coverage, theoretical_coverage
+    """
+    # load configuration file
+    if eiv:
+        with open(os.path.join('configurations',f'eiv_{data}.json'),'r') as\
+                conf_file:
+            conf_dict = json.load(conf_file)
+    else:
+        with open(os.path.join('configurations',f'noneiv_{data}.json'),'r') as\
+                conf_file:
+            conf_dict = json.load(conf_file)
+
+    long_dataname = conf_dict["long_dataname"]
+    short_dataname = conf_dict["short_dataname"]
+
+
+    load_data = importlib.import_module(f'EIVData.{long_dataname}').load_data
+
+    # switch to gpu, if possible
+    try:
+        gpu_number = conf_dict["gpu_number"]
+        device = torch.device(f'cuda:{gpu_number}')
+        try:
+            torch.tensor([0.0]).to(device)
+        except RuntimeError:
+            if torch.cuda.is_available():
+                print('Switched to GPU 0')
+                device = torch.device('cuda:0')
+            else:
+                print('No cuda available, using CPU')
+                device = torch.device('cpu')
+    except KeyError:
+        device = torch.device('cpu')
+
+
+    if use_ground_truth:
+        # test whether there is a ground truth
+        try:
+            train_data, _, _,_  \
+                    = load_data(seed=0, return_ground_truth=True)
+        except TypeError:
+        # if not, end function
+            return None,None
+    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()
+    output_dim = train_data[0][1].numel()
+
+    ## Create iterators for get_coverage_distribution
+    seed_list = range(conf_dict["seed_range"][0],
+            conf_dict["seed_range"][1])
+
+    # iterator for networks
+    def net_iterator(eiv=eiv, seed_list=seed_list):
+        """
+        Yields EiV models (if `eiv`) or
+        non-EiV models (if not `eiv`) for the seeds in
+        `seed_list` and `data`.
+        """
+        if eiv:
+            # load parameters
+            init_std_y = conf_dict["init_std_y_list"][0]
+            unscaled_reg = conf_dict["unscaled_reg"]
+            p = conf_dict["p"]
+            hidden_layers = conf_dict["hidden_layers"]
+            fixed_std_x = conf_dict["fixed_std_x"]
+            net = Networks.FNNEIV(p=p, init_std_y=init_std_y,
+                    h=[input_dim, *hidden_layers, output_dim],
+                    fixed_std_x=fixed_std_x).to(device)
+            for seed in seed_list:
+                # load network paramaters
+                saved_file = os.path.join('saved_networks',
+                        f'eiv_{short_dataname}'\
+                                f'_init_std_y_{init_std_y:.3f}'\
+                                f'_ureg_{unscaled_reg:.1f}'\
+                                f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                                f'_seed_{seed}.pkl')
+                train_and_store.open_stored_training(saved_file=saved_file,
+                        net=net, device=device)
+                yield net
+        else:
+            # load parameters
+            init_std_y = conf_dict["init_std_y_list"][0]
+            unscaled_reg = conf_dict["unscaled_reg"]
+            p = conf_dict["p"]
+            hidden_layers = conf_dict["hidden_layers"]
+            net = Networks.FNNBer(p=p, init_std_y=init_std_y,
+                    h=[input_dim, *hidden_layers, output_dim]).to(device)
+            for seed in seed_list:
+                saved_file = os.path.join('saved_networks',
+                            f'noneiv_{short_dataname}'\
+                                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                                    f'_p_{p:.2f}_seed_{seed}.pkl')
+                # load network paramaters
+                train_and_store.open_stored_training(saved_file=saved_file,
+                        net=net, device=device)
+                yield net
+
+    # iterator for dataloaders
+    def dataloader_iterator(seed_list=seed_list,use_ground_truth=use_ground_truth,
+            batch_size = 100):
+        """
+        Yields dataloaders for `data`, according to the seeds in `seed_list`.
+        If `use_ground_truth` the data is cut to contain noisy x and unnoisy y.
+        """
+        for seed in seed_list:
+            if not use_ground_truth:
+                _, test_data = load_data(seed=seed)
+                test_dataloader = DataLoader(test_data, 
+                        batch_size=batch_size,
+                        shuffle=True)
+                yield test_dataloader
+            else:
+                _, _, _, true_test =\
+                        load_data(seed=seed, return_ground_truth=True)
+                # take noisy x but unnoisy y
+                cut_true_test = VerticalCut(true_test,
+                        components_to_pick=[2,1])
+                test_dataloader = DataLoader(cut_true_test, 
+                        batch_size=batch_size,
+                        shuffle=True)
+                yield test_dataloader
+
+
+    # Compute coverages
+    numerical_coverage, theoretical_coverage = get_coverage_distribution(
+        net_iterator=net_iterator(eiv=eiv),
+        dataloader_iterator=dataloader_iterator(),
+        device=device,
+        number_of_draws=number_of_draws,
+        q_range=q_range,
+        noisy_y = not use_ground_truth)
+    return numerical_coverage, theoretical_coverage
+
+
+
+####
+
+def coverage_diagonal_plot(eiv_coverages, noneiv_coverages, color,
+        against_theoretical = False, label = '', mean_error=True):
+    """
+    Plot numerical coverages against q (used coverage value), if
+    against_theoretical is False, or the theoretical coverage, if
+    `against_theoretical` is True.
+    :param eiv_coverages: The output of `compute_coverages` with `eiv=True`
+    :param noneiv_coverages: The output of `compute_coverages` with `eiv=False`
+    :param color: String, denoting the color.
+    :param against_theoretical: Boolean, see above.
+    :param label: String. Will be included as label in the plot.
+    :param mean_error: Boolean. If True the standard deviation is divided
+    by the square root of the number of elements, to display the error 
+    of the mean (and not the dispersion).
+    """
+    eiv_numerical_coverage, eiv_theoretical_coverage = eiv_coverages
+    noneiv_numerical_coverage, noneiv_theoretical_coverage = noneiv_coverages
+    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 mean_error:
+        std_eiv_numerical_coverage /= np.sqrt(eiv_numerical_coverage.shape[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=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 mean_error:
+        std_noneiv_numerical_coverage /= \
+                np.sqrt(noneiv_numerical_coverage.shape[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')
+    # 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='',
+        mean_error=True):
+    """
+    Plot deviation of numerical coverages to q (used coverage value), if
+    against_theoretical is False, or to the theoretical coverage, if
+    `against_theoretical` is True. On the x axis q is plotted (in contrast to
+            coverage_diagonal_plot).
+    :param eiv_coverages: The output of `compute_coverages` with `eiv=True`
+    :param noneiv_coverages: The output of `compute_coverages` with `eiv=False`
+    :param color: String, denoting the color.
+    :param absolute_values:Boolean. If True (default) the absolute value of the
+    differences will be plotted.
+    :param against_theoretical: Boolean, see above.
+    :param label: String. Will be included as label in the plot.
+    :param mean_error: Boolean. If True the standard deviation is divided
+    by the square root of the number of elements, to display the error 
+    of the mean (and not the dispersion).
+    """
+    if not against_theoretical:
+        eiv_coverages = eiv_coverages[0]
+        noneiv_coverages = noneiv_coverages[0]
+        assert len(eiv_coverages.shape) == 2 and\
+                eiv_coverages.shape[0] == len(q_range)
+        assert len(eiv_coverages.shape) == 2 and\
+                eiv_coverages.shape[0] == len(q_range)
+        assert len(noneiv_coverages.shape) == 2
+        eiv_cov_res =  eiv_coverages - q_range.reshape((-1,1))
+        noneiv_cov_res = noneiv_coverages - q_range.reshape((-1,1))
+    else:
+        eiv_cov_res =  eiv_coverages[0] - eiv_coverages[1]
+        noneiv_cov_res =  noneiv_coverages[0] - noneiv_coverages[1]
+
+    if absolute_values:
+        eiv_cov_res = np.abs(eiv_cov_res)
+        noneiv_cov_res = np.abs(noneiv_cov_res)
+    mean_eiv_cov_res = np.mean(eiv_cov_res, axis=-1)
+    std_eiv_cov_res = np.std(eiv_cov_res, axis=-1)
+    if mean_error:
+        std_eiv_cov_res /= np.sqrt(eiv_cov_res.shape[1])
+    mean_noneiv_cov_res = np.mean(noneiv_cov_res, axis=-1)
+    std_noneiv_cov_res = np.std(noneiv_cov_res, axis=-1)
+    if mean_error:
+        std_noneiv_cov_res /= np.sqrt(noneiv_cov_res.shape[1])
+    plt.plot(q_range, mean_eiv_cov_res,
+            color=color, linestyle='solid', label=label)
+    plt.fill_between(q_range,
+            mean_eiv_cov_res - std_eiv_cov_res,
+            mean_eiv_cov_res + std_eiv_cov_res, 
+            color=color, alpha=0.5)
+    plt.plot(q_range, mean_noneiv_cov_res,
+            color=color, linestyle='dashed')
+    plt.fill_between(q_range,
+            mean_noneiv_cov_res - std_noneiv_cov_res,
+            mean_noneiv_cov_res + std_noneiv_cov_res,
+            color=color, alpha=0.3)
+
+
+def coverage_difference_plot(eiv_coverages, noneiv_coverages, color, label='',
+        mean_error=True):
+    """
+    Plot difference of numerical coverages from EiV to the one of non-EiV.
+    On the x axis q is plotted.
+    :param eiv_coverages: The output of `compute_coverages` with `eiv=True`
+    :param noneiv_coverages: The output of `compute_coverages` with `eiv=False`
+    :param color: String, denoting the color.
+    :param label: String. Will be included as label in the plot.
+    :param mean_error: Boolean. If True the standard deviation is divided
+    by the square root of the number of elements, to display the error 
+    of the mean (and not the dispersion).
+    """
+    cov_diff =  eiv_coverages[0] - noneiv_coverages[0]
+
+    mean_cov_diff = np.mean(cov_diff, axis=-1)
+    std_cov_diff = np.std(cov_diff, axis=-1)
+    if mean_error:
+        std_cov_diff /= np.sqrt(cov_diff.shape[1])
+    plt.plot(q_range, mean_cov_diff,
+            color=color, linestyle='solid', label=label)
+    plt.fill_between(q_range, mean_cov_diff - std_cov_diff,
+            mean_cov_diff + std_cov_diff, 
+            color=color, alpha=0.5)
+
+
+
+
+
+# create figures, together with title and axis labels
+plt.figure(1)
+plt.clf()
+plt.title('Coverage for datasets with ground truth')
+plt.xlabel('q')
+plt.ylabel('coverage')
+plt.figure(2)
+plt.clf()
+plt.title('Coverage of noisy labels (vs. theory)')
+plt.xlabel('th. coverage')
+plt.ylabel('coverage')
+plt.figure(3)
+plt.clf()
+plt.title('Coverage of noisy labels (vs. q)')
+plt.xlabel('q')
+plt.ylabel('coverage')
+plt.figure(4)
+plt.clf()
+plt.title('q Deviation of coverage of noisy labels')
+plt.xlabel('q')
+plt.ylabel('coverage')
+plt.figure(5)
+plt.clf()
+plt.title('q Deviation of coverage for datasets with ground truth')
+plt.xlabel('q')
+plt.ylabel('deviation cov from q')
+plt.figure(6)
+plt.clf()
+plt.title('Theory deviation  of noisy coverage')
+plt.xlabel('q')
+plt.ylabel('deviation cov from theor. cov')
+plt.figure(7)
+plt.clf()
+plt.title('Difference coverages')
+plt.xlabel('q')
+plt.ylabel('EiV cov - nonEiV-cov')
+
+# datasets to plot and their coloring
+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]:
+    for data, color in zip(datasets, colors):
+        # compute coverages
+        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)
+        # If no result was produced (no ground truth although required), skip
+        if eiv_coverages[0] is None and eiv_coverages[1] is None:
+            assert noneiv_coverages[0] is None and noneiv_coverages[1] is None
+            break
+        
+        # create plots
+        if use_ground_truth:
+            plt.figure(1)
+            coverage_diagonal_plot(eiv_coverages, noneiv_coverages, 
+                    color=color, against_theoretical=False, label=data)
+            plt.figure(5)
+            coverage_residual_plot(eiv_coverages, noneiv_coverages, 
+                    color=color, against_theoretical=False, label=data)
+        else:
+            plt.figure(2)
+            coverage_diagonal_plot(eiv_coverages, noneiv_coverages, 
+                    color=color, against_theoretical=True, label=data)
+            plt.figure(3)
+            coverage_diagonal_plot(eiv_coverages, noneiv_coverages, 
+                    color=color, against_theoretical=False, label=data)
+            plt.figure(4)
+            coverage_residual_plot(eiv_coverages, noneiv_coverages, 
+                    color=color, against_theoretical=False, label=data)
+            plt.figure(6)
+            coverage_residual_plot(eiv_coverages, noneiv_coverages, 
+                    color=color, against_theoretical=True, label=data)
+            plt.figure(7)
+            coverage_difference_plot(eiv_coverages, noneiv_coverages, 
+                    color=color, label=data)
+
+# add diagonals, where meaningful
+plt.figure(1)
+x_diag = np.linspace(0.0, 1.0)
+plt.plot(x_diag, x_diag, color='k', linestyle='dotted' )
+plt.figure(2)
+x_diag = np.linspace(0.0, 1.0)
+plt.plot(x_diag, x_diag, color='k', linestyle='dotted' )
+plt.figure(3)
+x_diag = np.linspace(0.0, 1.0)
+plt.plot(x_diag, x_diag, color='k', linestyle='dotted' )
+
+
+# add legends
+for fig_nr in range(1,8):
+    plt.figure(fig_nr)
+    plt.legend()
+
+plt.figure(1)
+plt.savefig('results/figures/summary_coverage_ground_truth.pdf')
+plt.figure(2)
+plt.savefig('results/figures/summary_coverage_noisy_vs_th.pdf')
+plt.figure(3)
+plt.savefig('results/figures/summary_coverage_noisy_vs_q.pdf')
+plt.figure(4)
+plt.savefig('results/figures/summary_q_deviation_coverage_noisy.pdf')
+plt.figure(5)
+plt.savefig('results/figures/summary_q_deviation_coverage_ground_truth.pdf')
+plt.figure(6)
+plt.savefig('results/figures/summary_th_deviation_coverage_ground_truth.pdf')
+plt.figure(7)
+plt.savefig('results/figures/summary_difference_coverage.pdf')
+
+# show plots
+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