diff --git a/EIVPackage/EIVArchitectures/Networks.py b/EIVPackage/EIVArchitectures/Networks.py
index 2d8186563e6d669ba77603e77e507c411260910f..09d75b877ba2617af15447e9bacd22b8581de35d 100644
--- a/EIVPackage/EIVArchitectures/Networks.py
+++ b/EIVPackage/EIVArchitectures/Networks.py
@@ -14,7 +14,7 @@ class FNNEIV(nn.Module):
     """
     A fully connected net with Error-in-Variables input and Bernoulli dropout
     layers. 
-    :param p: dropout rate, defaults to 0.5
+    :param p: dropout rate, defaults to 0.2
     :param init_std_y: Initial estimated standard deviation for y.
     :param precision_prior_zeta: precision of the prior for zeta.
     Defaults to 0.0 (=improper prior)
@@ -23,7 +23,10 @@ class FNNEIV(nn.Module):
     `fixed_std_x` is different from `None`.
     :param h: A list specifying the number of neurons in each layer.
     :param fixed_std_x: If given, this value will be the output of the method
-    `get_std_x()`.
+    `get_std_x()`, no matter what the deming factor.
+    :param repetition: Positive integer, the default value for repeating input,
+    defaults to 1.  For a single call this can also be specified in the forward
+    method.
     **Note**: 
     - To change the deming factor afterwards, use the method `change_deming`
     - To change fixed_std_x afterwards, use the method `change_fixed_std_x`
@@ -31,14 +34,14 @@ class FNNEIV(nn.Module):
     LeakyReLUSlope = 1e-2
     def __init__(self, p = 0.2, init_std_y=1.0, precision_prior_zeta=0.0, 
             deming=1.0, h=[10, 1024,1024,1024,1024, 1], 
-            fixed_std_x = None):
+            fixed_std_x = None, repetition = 1):
         super().__init__()
         # part before Bernoulli dropout
         self.init_std_y = init_std_y
         InverseSoftplus = lambda sigma: torch.log(torch.exp(sigma) - 1 )
-        self.std_y_par = nn.Parameter(
+        self.std_y_par = nn.parameter.Parameter(
                 InverseSoftplus(torch.tensor([init_std_y])))
-        self._repetition = 1
+        self._repetition = repetition
         self.main = nn.Sequential(
                 EIVInput(precision_prior_zeta=precision_prior_zeta, 
                     external_std_x=self.get_std_x),
@@ -57,6 +60,9 @@ class FNNEIV(nn.Module):
                 nn.Linear(h[4], h[5]))
         self.p = p
         self._deming = deming
+        if fixed_std_x is not None:
+            if type(fixed_std_x) is not torch.tensor:
+                fixed_std_x = torch.tensor(fixed_std_x)
         self._fixed_std_x = fixed_std_x
         # needed for switch_noise_off()
         self.noise_is_on = True
@@ -76,6 +82,9 @@ class FNNEIV(nn.Module):
         :param fixed_std_x: A positive float
         """
         print('Updating fixed_std_x from %.3f to %.3f' % (self._fixed_std_x, fixed_std_x))
+        if fixed_std_x is not None:
+            if type(fixed_std_x) is not torch.tensor:
+                fixed_std_x = torch.tensor(fixed_std_x)
         self._fixed_std_x = fixed_std_x
 
     def noise_off(self):
@@ -95,7 +104,7 @@ class FNNEIV(nn.Module):
             else:
                 return self._fixed_std_x
         else:
-            return 0.0
+            return torch.tensor(0.0, dtype=torch.float32)
 
     def get_std_y(self):
         return nn.Softplus()(self.std_y_par)
@@ -104,10 +113,11 @@ class FNNEIV(nn.Module):
         return self._repetition
 
     def forward(self, x, repetition=1):
+        old_repetition = self._repetition
         self._repetition = repetition
         mu = self.main(x)
         sigma = self.sigma(mu)
-        self._repetition = 1
+        self._repetition = old_repetition
         return mu, sigma
 
     def regularizer(self, x, lamb):
@@ -148,8 +158,10 @@ class FNNEIV(nn.Module):
                 regularization += self.main[0].neg_x_evidence(x)
         return regularization 
 
-
-    def predict(self, x, number_of_draws=100, remove_graph=True
+    # TODO: repetition
+    # + Test via breakpoints
+    def predict(self, x, number_of_draws=[100,5], number_of_parameter_chunks = None,
+            remove_graph=True
             , take_average_of_prediction=True):
         """
         Average over `number_of_draws` forward passes. If 
@@ -158,19 +170,46 @@ class FNNEIV(nn.Module):
         **Note**: This method does neither touch the input noise nor Dropout.
         The corresponding setting is left to the user!
         :param x: A torch.tensor, the input
-        :param number_of_draws: Number of draws to obtain from x
+        :param number_of_draws: An integer or a list. If an integer
+        `number_of_draws`, will be converted internally to
+        `[number_of_draws,1]`.Numbers of draws to obtain from x via parameter
+        sampling (first element) and noise input sampling (second element).
+        :param number_of_parameter_chunks: An integer or None (default). If
+        None, the second element of `number_of_draws` will be taken (and will
+        thus be identical to 1 if `number_of_draws` is an integer, see above).
+        Samples in the parameter space will be divided into
+        `number_of_parameter_chunks` chunks when collected. Can be used to
+        reduced the memory usage. 
         :param remove_graph: If True (default) the output will 
         be detached to save memory
         :param take_average_of_prediction: If False, no averaging will be 
         applied to the prediction and the second dimension of the first output
         will count the number_of_draws.
         """
-        x, = repeat_tensors(x, number_of_draws=number_of_draws)
-        pred, sigma = self.forward(x)
-        if remove_graph:
-            pred, sigma = pred.detach(), sigma.detach()
-        pred, sigma = reshape_to_chunks(pred, sigma, 
-                number_of_draws=number_of_draws)
+        # x, = repeat_tensors(x, number_of_draws=number_of_draws)
+        if type(number_of_draws) is int:
+            number_of_draws = [number_of_draws, 1]
+        if number_of_parameter_chunks is None:
+            number_of_parameter_chunks = number_of_draws[1]
+        chunk_size = number_of_draws[0] // number_of_parameter_chunks
+        pred_collection, sigma_collection = [], []
+        remaining_draws = number_of_draws[0]
+        while remaining_draws > 0:
+            if remaining_draws < chunk_size:
+                parameter_sample_size = remaining_draws
+            else:
+                parameter_sample_size = chunk_size
+            repeated_x, = repeat_tensors(x, number_of_draws=parameter_sample_size * number_of_draws[1])
+            pred, sigma = self.forward(repeated_x, repetition=number_of_draws[1])
+            if remove_graph:
+                pred, sigma = pred.detach(), sigma.detach()
+            pred, sigma = reshape_to_chunks(pred, sigma, 
+                    number_of_draws=parameter_sample_size * number_of_draws[1])
+            pred_collection.append(pred)
+            sigma_collection.append(sigma)
+            remaining_draws -= parameter_sample_size
+        pred = torch.cat(pred_collection, dim=1)
+        sigma = torch.cat(sigma_collection, dim=1)
         # reduce along the draws (actually redundant for sigma)
         if take_average_of_prediction:
             pred, sigma = torch.mean(pred, dim=1), torch.mean(sigma, dim=1)
@@ -178,6 +217,75 @@ class FNNEIV(nn.Module):
             sigma = torch.mean(sigma, dim=1)
         return pred, sigma
 
+    # TODO: repetition
+    # + Implement
+    # + Test via breakpoints
+    def predictive_logdensity(self, x, y, number_of_draws=[100, 5], number_of_parameter_chunks = None, remove_graph=True,
+            average_batch_dimension=True, scale_labels=None,
+            decouple_dimensions=False):
+        """
+        Computes the logarithm of the predictive density evaluated at `y`. If
+        `average_batch_dimension` is `True` these values will be averaged over
+        the batch dimension.
+        :param x: A torch.tensor, the input
+        :param y: A torch.tensor, labels on which to evaluate the density
+        :param number_of_draws: An integer or a list. If an integer
+        `number_of_draws`, will be converted internally to
+        `[number_of_draws,1]`.Numbers of draws to obtain from x via parameter
+        sampling (first element) and noise input sampling (second element).
+        :param number_of_parameter_chunks: An integer or None (default). If
+        None, the second element of `number_of_draws` will be taken (and will
+        thus be identical to 1 if `number_of_draws` is an integer, see above).
+        Samples in the parameter space will be divided into
+        `number_of_parameter_chunks` chunks when collected. Can be used to
+        reduced the memory usage. 
+        :param remove_graph: If True (default) the output will 
+        be detached to save memory
+        :param average_batch_dimension: Boolean. If True (default) the values
+        will be averaged over the batch dimension. If False, the batch
+        dimension will be left untouched and all values will be returned.
+        """
+        out, sigmas = self.predict(x, number_of_draws=number_of_draws,
+                number_of_parameter_chunks=number_of_parameter_chunks,
+                remove_graph=remove_graph,
+                take_average_of_prediction=False)
+        # Add "repetition" dimension to y and out
+        y = y[:,None,...]
+        sigmas = sigmas[:,None,...]
+        if len(y.shape) <= 2:
+            # add an output axis if necessary
+            y = y[...,None]
+            sigmas = sigmas[...,None]
+        # squeeze last dimensions into one
+        y = y.view((*y.shape[:2], -1))
+        sigmas = sigmas.view((*sigmas.shape[:2], -1))
+        out = out.view((*out.shape[:2], -1))
+        # check if dimensions consistent
+        assert y.shape == sigmas.shape
+        assert y.shape[0] == out.shape[0]
+        assert y.shape[2] == out.shape[2]
+        if scale_labels is not None:
+            extended_scale_labels = scale_labels.flatten()[None,None,:]
+            out = out * extended_scale_labels
+            y = y * extended_scale_labels
+            sigmas = sigmas * extended_scale_labels
+        # exponential argument for density
+        if not decouple_dimensions:
+            exp_arg =  torch.sum(-1/(2*sigmas**2) * (y-out)**2-\
+                        1/2 * torch.log(2 * torch.pi * sigmas**2), dim=2)
+        else:
+            exp_arg =  -1/(2*sigmas**2) * (y-out)**2-\
+                            1/2 * torch.log(2 * torch.pi * sigmas**2)
+        # average over parameter values
+        predictive_log_density_values = \
+                torch.logsumexp(input=exp_arg, dim=1)\
+                    - torch.log(torch.prod(torch.tensor(number_of_draws))) 
+        if average_batch_dimension:
+            return torch.mean(predictive_log_density_values, dim=0)
+        else:
+            return predictive_log_density_values
+
+
 class FNNBer(nn.Module):
     """
     A fully connected net Bernoulli dropout layers.
@@ -191,7 +299,7 @@ class FNNBer(nn.Module):
         # part before Bernoulli dropout
         self.init_std_y = init_std_y
         InverseSoftplus = lambda sigma: torch.log(torch.exp(sigma) - 1 )
-        self.std_y_par = nn.Parameter(
+        self.std_y_par = nn.parameter.Parameter(
                 InverseSoftplus(torch.tensor([init_std_y])))
         self.main = nn.Sequential(
                 nn.Linear(h[0], h[1]),
@@ -263,8 +371,9 @@ class FNNBer(nn.Module):
         :param remove_graph: If True (default) the output will 
         be detached to save memory
         :param take_average_of_prediction: If False, no averaging will be 
-        applied to the prediction and the second dimension of the first output
+        applied to the prediction and the second dimension of the first output 
         will count the number_of_draws.
+        :returns: predictions, sigmas
         """
         x, = repeat_tensors(x, number_of_draws=number_of_draws)
         pred, sigma = self.forward(x)
diff --git a/EIVPackage/EIVData/california_housing.py b/EIVPackage/EIVData/california_housing.py
index d1f442413324186ab6e51bfa59049e1cb5209e96..cf755199ff8c1a76d03fc4c17c4866a61084dbb2 100644
--- a/EIVPackage/EIVData/california_housing.py
+++ b/EIVPackage/EIVData/california_housing.py
@@ -1,3 +1,4 @@
+import torch
 from EIVData.csv_dataset import CSVData
 from torch.utils.data import random_split
 
@@ -18,7 +19,9 @@ def load_data(seed=0, splitting_part=0.8, normalize=True):
     dataset_len = len(california_dataset)
     train_len = int(dataset_len*splitting_part)
     test_len = dataset_len - train_len
-    california_trainset, california_testset = random_split(california_dataset , lengths=[train_len, test_len])
+    california_trainset, california_testset = random_split(california_dataset,
+            lengths=[train_len, test_len],
+            generator=torch.Generator().manual_seed(seed))
     return california_trainset, california_testset
 
 
diff --git a/EIVPackage/EIVData/concrete_strength.py b/EIVPackage/EIVData/concrete_strength.py
index 472548c9958758d83cead8f6930a918d0c241e16..d6d13c4446dc8f5fd3b52690d7382d9dc8ba3b25 100644
--- a/EIVPackage/EIVData/concrete_strength.py
+++ b/EIVPackage/EIVData/concrete_strength.py
@@ -1,3 +1,4 @@
+import torch
 from EIVData.csv_dataset import CSVData
 from torch.utils.data import random_split
 
@@ -11,12 +12,14 @@ def load_data(seed=0, splitting_part=0.8, normalize=True):
     :normalize: Whether to normalize the data, defaults to True.
     :returns: concrete_trainset, concrete_testset
     """
-    concrete_dataset = CSVData('~/SharedData/AI/datasets/concrete_compression_strength/compresive_strength_concrete.csv',
+    concrete_dataset = CSVData('~/SharedData/AI/datasets/concrete_compression_strength/compressive_strength_concrete.csv',
             class_name='Concrete compressive strength(MPa, megapascals) ',
             shuffle_seed=seed,
             normalize=normalize)
     dataset_len = len(concrete_dataset)
     train_len = int(dataset_len*splitting_part)
     test_len = dataset_len - train_len
-    concrete_trainset, concrete_testset = random_split(concrete_dataset , lengths=[train_len, test_len])
+    concrete_trainset, concrete_testset = random_split(concrete_dataset,
+            lengths=[train_len, test_len],
+            generator=torch.Generator().manual_seed(seed))
     return concrete_trainset, concrete_testset 
diff --git a/EIVPackage/EIVData/energy_efficiency.py b/EIVPackage/EIVData/energy_efficiency.py
index b421b9c03cad95ba5271372f7c2d76d13358617f..ddf63170c2ddf2beb1d2aa96d85059f50798c27f 100644
--- a/EIVPackage/EIVData/energy_efficiency.py
+++ b/EIVPackage/EIVData/energy_efficiency.py
@@ -1,3 +1,4 @@
+import torch
 from EIVData.csv_dataset import CSVData
 from torch.utils.data import random_split
 
@@ -18,5 +19,7 @@ def load_data(seed=0, splitting_part=0.8, normalize=True):
     dataset_len = len(energy_dataset)
     train_len = int(dataset_len*splitting_part)
     test_len = dataset_len - train_len
-    energy_trainset, energy_testset = random_split(energy_dataset , lengths=[train_len, test_len])
+    energy_trainset, energy_testset = random_split(energy_dataset,
+            lengths=[train_len, test_len],
+            generator=torch.Generator().manual_seed(seed))
     return energy_trainset, energy_testset
diff --git a/EIVPackage/EIVData/kin8nm.py b/EIVPackage/EIVData/kin8nm.py
index 7cb7ad3b3b3d782f8fd8e15822d1dd5e3752e4a9..774f1f46c9549ba92754b5dd71870ca70796185d 100644
--- a/EIVPackage/EIVData/kin8nm.py
+++ b/EIVPackage/EIVData/kin8nm.py
@@ -1,3 +1,4 @@
+import torch
 from EIVData.csv_dataset import CSVData
 from torch.utils.data import random_split
 
@@ -18,7 +19,9 @@ def load_data(seed=0, splitting_part=0.8, normalize=True):
     dataset_len = len(kin8nm_dataset)
     train_len = int(dataset_len*splitting_part)
     test_len = dataset_len - train_len
-    kin8nm_trainset, kin8nm_testset = random_split(kin8nm_dataset , lengths=[train_len, test_len])
+    kin8nm_trainset, kin8nm_testset = random_split(kin8nm_dataset,
+            lengths=[train_len, test_len],
+            generator=torch.Generator().manual_seed(seed))
     return kin8nm_trainset, kin8nm_testset
 
 
diff --git a/EIVPackage/EIVData/million_song.py b/EIVPackage/EIVData/million_song.py
index 1dd545ac7297157b2596e4613278b615bdc8869a..ae6d18a446973a3c33f5a26289bf1d044827f340 100644
--- a/EIVPackage/EIVData/million_song.py
+++ b/EIVPackage/EIVData/million_song.py
@@ -1,3 +1,4 @@
+import torch
 from EIVData.csv_dataset import CSVData
 from torch.utils.data import random_split
 
@@ -21,5 +22,6 @@ def load_data(seed=0, splitting_part=0.8, normalize=True):
     train_len = int(dataset_len*splitting_part)
     test_len = dataset_len - train_len
     msd_trainset, msd_testset = random_split(msd_dataset,
-            lengths=[train_len, test_len])
+            lengths=[train_len, test_len],
+            generator=torch.Generator().manual_seed(seed))
     return msd_trainset, msd_testset 
diff --git a/EIVPackage/EIVData/naval_propulsion.py b/EIVPackage/EIVData/naval_propulsion.py
index a3e95962324bf5698a88325494741d44d172d7f0..8749899f76108888085e1361d3da7940d7cd21d9 100644
--- a/EIVPackage/EIVData/naval_propulsion.py
+++ b/EIVPackage/EIVData/naval_propulsion.py
@@ -1,3 +1,4 @@
+import torch
 from EIVData.csv_dataset import CSVData
 from torch.utils.data import random_split
 
@@ -20,6 +21,8 @@ def load_data(seed=0, splitting_part=0.8, normalize=True):
     dataset_len = len(naval_dataset)
     train_len = int(dataset_len*splitting_part)
     test_len = dataset_len - train_len
-    naval_trainset, naval_testset = random_split(naval_dataset , lengths=[train_len, test_len])
+    naval_trainset, naval_testset = random_split(naval_dataset,
+            lengths=[train_len, test_len],
+            generator=torch.Generator().manual_seed(seed))
     return naval_trainset, naval_testset
 
diff --git a/EIVPackage/EIVData/power_plant.py b/EIVPackage/EIVData/power_plant.py
index e33f7420dd716860c63ee57eb7b7e18e0b2f371a..6f40c65d0ac2b201f2549a347f9823cba9a10035 100644
--- a/EIVPackage/EIVData/power_plant.py
+++ b/EIVPackage/EIVData/power_plant.py
@@ -1,3 +1,4 @@
+import torch
 from EIVData.csv_dataset import CSVData
 from torch.utils.data import random_split
 
@@ -19,6 +20,8 @@ def load_data(seed=0, splitting_part=0.8, normalize=True):
     dataset_len = len(naval_dataset)
     train_len = int(dataset_len*splitting_part)
     test_len = dataset_len - train_len
-    naval_trainset, naval_testset = random_split(naval_dataset , lengths=[train_len, test_len])
+    naval_trainset, naval_testset = random_split(naval_dataset,
+            lengths=[train_len, test_len],
+            generator=torch.Generator().manual_seed(seed))
     return naval_trainset, naval_testset
 
diff --git a/EIVPackage/EIVData/protein_structure.py b/EIVPackage/EIVData/protein_structure.py
index 50e34e1716405fde103527562bae4da61a2b8c2a..83050871a5d2d4c1411d4ae4053a877a5cad153e 100644
--- a/EIVPackage/EIVData/protein_structure.py
+++ b/EIVPackage/EIVData/protein_structure.py
@@ -1,3 +1,4 @@
+import torch
 from EIVData.csv_dataset import CSVData
 from torch.utils.data import random_split
 
@@ -20,5 +21,6 @@ def load_data(seed=0, splitting_part=0.8, normalize=True):
     train_len = int(dataset_len*splitting_part)
     test_len = dataset_len - train_len
     protein_trainset, protein_testset = random_split(protein_dataset,
-            lengths=[train_len, test_len])
+            lengths=[train_len, test_len],
+            generator=torch.Generator().manual_seed(seed))
     return protein_trainset, protein_testset 
diff --git a/EIVPackage/EIVData/wine_quality.py b/EIVPackage/EIVData/wine_quality.py
index eb4216aa2a25a6d16c72345496a81164f8c08920..78c839a2fc23e6ad8a315f985271303a6cd4ae62 100644
--- a/EIVPackage/EIVData/wine_quality.py
+++ b/EIVPackage/EIVData/wine_quality.py
@@ -1,3 +1,4 @@
+import torch
 from EIVData.csv_dataset import CSVData
 from torch.utils.data import random_split
 
@@ -20,5 +21,6 @@ def load_data(seed=0, splitting_part=0.8, normalize=True):
     train_len = int(dataset_len*splitting_part)
     test_len = dataset_len - train_len
     wine_trainset, wine_testset = random_split(wine_dataset,
-            lengths=[train_len, test_len])
+            lengths=[train_len, test_len],
+            generator=torch.Generator().manual_seed(seed))
     return wine_trainset, wine_testset 
diff --git a/EIVPackage/EIVData/yacht_hydrodynamics.py b/EIVPackage/EIVData/yacht_hydrodynamics.py
index f3cf3f8e126e51feae21387ba8bb1c8c6a718b9e..505c841ac8e369bbac1c547e1d351b92641e1f49 100644
--- a/EIVPackage/EIVData/yacht_hydrodynamics.py
+++ b/EIVPackage/EIVData/yacht_hydrodynamics.py
@@ -1,3 +1,4 @@
+import torch
 from EIVData.csv_dataset import CSVData
 from torch.utils.data import random_split
 
@@ -21,7 +22,8 @@ def load_data(seed=0, splitting_part=0.8, normalize=True):
     train_len = int(dataset_len*splitting_part)
     test_len = dataset_len - train_len
     yacht_trainset, yacht_testset = random_split(yacht_dataset,
-            lengths=[train_len, test_len])
+            lengths=[train_len, test_len],
+            generator=torch.Generator().manual_seed(seed))
     return yacht_trainset, yacht_testset
 
 
diff --git a/EIVPackage/EIVTrainingRoutines/loss_functions.py b/EIVPackage/EIVTrainingRoutines/loss_functions.py
index 37b76fe46f0518708a8f6fc45c38d2784d604db1..88e924df76c1b2bd00e4c4abcdf18089f699539e 100644
--- a/EIVPackage/EIVTrainingRoutines/loss_functions.py
+++ b/EIVPackage/EIVTrainingRoutines/loss_functions.py
@@ -17,6 +17,7 @@ def nll_reg_loss(net, x, y, reg):
     :param reg: A non-negative float, the regularization.
     """
     out, std_y = net(x)
+    # Add label dimension to y if missing
     if len(y.shape) <= 1:
         y = y.view((-1,1))
     assert out.shape == y.shape
@@ -26,13 +27,11 @@ def nll_reg_loss(net, x, y, reg):
     return neg_log_likelihood + regularization
 
 
-def nll_eiv_no_jensen(net, x, y, reg, number_of_draws=5):
+def nll_eiv(net, x, y, reg, number_of_draws=5):
     """
     negative log likelihood criterion for an Error in variables model (EIV) 
     where `torch.logsumexp` is applied to partitions of size `number_of_draws`
     of `mu` and `sigma` in the batch dimension (that is the first one). 
-    **Note**: This function is supposed to be used in combination 
-    of `repeat_tensors` with the same argument `number_of_draws`.
     *Note that `reg` will not be divided by the data size (and by 2), 
      this should be done beforehand.*
     :param mu: predicted mu
@@ -40,13 +39,16 @@ def nll_eiv_no_jensen(net, x, y, reg, number_of_draws=5):
     :param y: ground truth
     :number_of_draws: Integer, supposed to be larger than 2
     """
+    # Add label dimension to y if missing
+    if len(y.shape) <= 1:
+        y = y.view((-1,1))
     regularization = net.regularizer(x, lamb=reg)
     # repeat_tensors
     x, y = repeat_tensors(x, y, number_of_draws=number_of_draws)
     pred, sigma = net(x, repetition=number_of_draws) 
     # split into chunks of size number_of_draws along batch dimension
-    pred, sigma, y = reshape_to_chunks(pred, sigma, 
-            y, number_of_draws=number_of_draws)
+    pred, sigma, y = reshape_to_chunks(pred, sigma, y, number_of_draws=number_of_draws)
+    assert pred.shape == y.shape
     # apply logsumexp to chunks and average the results
     nll = -1 * (torch.logsumexp(-1 * sigma.log()
         -((y-pred)**2)/(2*sigma**2), dim=1)
diff --git a/Experiments/evaluate_california.py b/Experiments/evaluate_california.py
deleted file mode 100644
index be631efaf4ca6cd6851eebeac5a5c5dcc7d13e8d..0000000000000000000000000000000000000000
--- a/Experiments/evaluate_california.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import os
-import numpy as np
-import torch
-import torch.backends.cudnn
-from torch.utils.data import DataLoader
-from torch.utils.tensorboard.writer import SummaryWriter
-
-from EIVArchitectures import Networks, initialize_weights
-from EIVData.california_housing import load_data
-from EIVTrainingRoutines import train_and_store, loss_functions
-
-from train_noneiv_california import p, init_std_y_list, seed_list, unscaled_reg, hidden_layers
-
-
-train_data, test_data = load_data()
-test_dataloader = DataLoader(test_data, batch_size=int(np.max((len(test_data), 800))))
-
-seed = seed_list[0]
-init_std_y = init_std_y_list[0]
-saved_file = os.path.join('saved_networks',
-            f'noneiv_california'\
-                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
-                    f'_p_{p:.2f}_seed_{seed}.pkl')
-
-input_dim = train_data[0][0].numel()
-output_dim = train_data[0][1].numel()
-net = Networks.FNNBer(p=p, init_std_y=init_std_y,
-        h=[input_dim, *hidden_layers, output_dim])
-train_and_store.open_stored_training(saved_file=saved_file,
-        net=net)
-
-
-# RMSE
-x,y = next(iter(test_dataloader))
-out = net(x)[0]
-if len(y.shape) <=1:
-    y = y.view((-1,1))
-assert y.shape == out.shape
-res = y-out
-scale = train_data.dataset.std_labels
-scaled_res = res * scale.view((1,-1))
-scaled_res = scaled_res.detach().cpu().numpy().flatten()
-rmse = np.sqrt(np.mean(scaled_res**2)) 
-print(f'RMSE {rmse:.3f}')
-
-
-# NLL
-x,y = next(iter(test_dataloader))
-training_state = net.training
-net.train()
-logdens = net.predictive_logdensity(x, y, number_of_draws=100,
-        decouple_dimensions=True,
-        scale_labels=train_data.dataset.std_labels.view((-1,))).mean()
-if training_state:
-    net.train()
-else:
-    net.eval()
-print(f'Dropout predictive {logdens:.3f}')
diff --git a/Experiments/evaluate_energy.py b/Experiments/evaluate_energy.py
deleted file mode 100644
index e9d74e3d6f9351d54221a41ef51c7b8d4f4f318c..0000000000000000000000000000000000000000
--- a/Experiments/evaluate_energy.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import os
-import numpy as np
-import torch
-import torch.backends.cudnn
-from torch.utils.data import DataLoader
-from torch.utils.tensorboard.writer import SummaryWriter
-
-from EIVArchitectures import Networks, initialize_weights
-from EIVData.energy_efficiency import load_data
-from EIVTrainingRoutines import train_and_store, loss_functions
-
-from train_noneiv_energy import p, init_std_y_list, seed_list, unscaled_reg, hidden_layers
-
-
-train_data, test_data = load_data()
-test_dataloader = DataLoader(test_data, batch_size=int(np.max((len(test_data), 800))))
-
-seed = seed_list[0]
-init_std_y = init_std_y_list[0]
-saved_file = os.path.join('saved_networks',
-            f'noneiv_energy'\
-                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
-                    f'_p_{p:.2f}_seed_{seed}.pkl')
-
-input_dim = train_data[0][0].numel()
-output_dim = train_data[0][1].numel()
-net = Networks.FNNBer(p=p, init_std_y=init_std_y,
-        h=[input_dim, *hidden_layers, output_dim])
-train_and_store.open_stored_training(saved_file=saved_file,
-        net=net)
-
-
-# RMSE
-x,y = next(iter(test_dataloader))
-out = net(x)[0]
-if len(y.shape) <=1:
-    y = y.view((-1,1))
-assert y.shape == out.shape
-res = y-out
-scale = train_data.dataset.std_labels
-scaled_res = res * scale.view((1,-1))
-scaled_res = scaled_res.detach().cpu().numpy().flatten()
-rmse = np.sqrt(np.mean(scaled_res**2)) 
-print(f'RMSE {rmse:.3f}')
-
-
-# NLL
-x,y = next(iter(test_dataloader))
-training_state = net.training
-net.train()
-logdens = net.predictive_logdensity(x, y, number_of_draws=100,
-        decouple_dimensions=True,
-        scale_labels=train_data.dataset.std_labels.view((-1,))).mean()
-if training_state:
-    net.train()
-else:
-    net.eval()
-print(f'Dropout predictive {logdens:.3f}')
diff --git a/Experiments/evaluate_kin8nm.py b/Experiments/evaluate_kin8nm.py
deleted file mode 100644
index a9d8ae6ab77700b1f32b127c89c5a2f9cf7c0ae2..0000000000000000000000000000000000000000
--- a/Experiments/evaluate_kin8nm.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import os
-import numpy as np
-import torch
-import torch.backends.cudnn
-from torch.utils.data import DataLoader
-from torch.utils.tensorboard.writer import SummaryWriter
-
-from EIVArchitectures import Networks, initialize_weights
-from EIVData.kin8nm import load_data
-from EIVTrainingRoutines import train_and_store, loss_functions
-
-from train_noneiv_kin8nm import p, init_std_y_list, seed_list, unscaled_reg, hidden_layers
-
-
-train_data, test_data = load_data()
-test_dataloader = DataLoader(test_data, batch_size=int(np.max((len(test_data), 800))))
-
-seed = seed_list[0]
-init_std_y = init_std_y_list[0]
-saved_file = os.path.join('saved_networks',
-            f'noneiv_kin8nm'\
-                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
-                    f'_p_{p:.2f}_seed_{seed}.pkl')
-
-input_dim = train_data[0][0].numel()
-output_dim = train_data[0][1].numel()
-net = Networks.FNNBer(p=p, init_std_y=init_std_y,
-        h=[input_dim, *hidden_layers, output_dim])
-train_and_store.open_stored_training(saved_file=saved_file,
-        net=net)
-
-
-# RMSE
-x,y = next(iter(test_dataloader))
-out = net(x)[0]
-if len(y.shape) <=1:
-    y = y.view((-1,1))
-assert y.shape == out.shape
-res = y-out
-scale = train_data.dataset.std_labels
-scaled_res = res * scale.view((1,-1))
-scaled_res = scaled_res.detach().cpu().numpy().flatten()
-rmse = np.sqrt(np.mean(scaled_res**2)) 
-print(f'RMSE {rmse:.3f}')
-
-
-# NLL
-x,y = next(iter(test_dataloader))
-training_state = net.training
-net.train()
-logdens = net.predictive_logdensity(x, y, number_of_draws=100,
-        decouple_dimensions=True,
-        scale_labels=train_data.dataset.std_labels.view((-1,))).mean()
-if training_state:
-    net.train()
-else:
-    net.eval()
-print(f'Dropout predictive {logdens:.3f}')
diff --git a/Experiments/evaluate_msd.py b/Experiments/evaluate_msd.py
deleted file mode 100644
index 041b7556da3eb8d7f206284f20ec0baf4eec2167..0000000000000000000000000000000000000000
--- a/Experiments/evaluate_msd.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import os
-import numpy as np
-import torch
-import torch.backends.cudnn
-from torch.utils.data import DataLoader
-from torch.utils.tensorboard.writer import SummaryWriter
-
-from EIVArchitectures import Networks, initialize_weights
-from EIVData.million_song import load_data
-from EIVTrainingRoutines import train_and_store, loss_functions
-
-from train_noneiv_msd import p, init_std_y_list, seed_list, unscaled_reg, hidden_layers
-
-
-train_data, test_data = load_data()
-test_dataloader = DataLoader(test_data, batch_size=int(np.max((len(test_data), 800))))
-
-seed = seed_list[0]
-init_std_y = init_std_y_list[0]
-saved_file = os.path.join('saved_networks',
-            f'noneiv_msd'\
-                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
-                    f'_p_{p:.2f}_seed_{seed}.pkl')
-
-input_dim = train_data[0][0].numel()
-output_dim = train_data[0][1].numel()
-net = Networks.FNNBer(p=p, init_std_y=init_std_y,
-        h=[input_dim, *hidden_layers, output_dim])
-train_and_store.open_stored_training(saved_file=saved_file,
-        net=net)
-
-
-# RMSE
-x,y = next(iter(test_dataloader))
-out = net(x)[0]
-if len(y.shape) <=1:
-    y = y.view((-1,1))
-assert y.shape == out.shape
-res = y-out
-scale = train_data.dataset.std_labels
-scaled_res = res * scale.view((1,-1))
-scaled_res = scaled_res.detach().cpu().numpy().flatten()
-rmse = np.sqrt(np.mean(scaled_res**2)) 
-print(f'RMSE {rmse:.3f}')
-
-
-# NLL
-x,y = next(iter(test_dataloader))
-training_state = net.training
-net.train()
-logdens = net.predictive_logdensity(x, y, number_of_draws=100,
-        decouple_dimensions=True,
-        scale_labels=train_data.dataset.std_labels.view((-1,))).mean()
-if training_state:
-    net.train()
-else:
-    net.eval()
-print(f'Dropout predictive {logdens:.3f}')
diff --git a/Experiments/evaluate_naval.py b/Experiments/evaluate_naval.py
deleted file mode 100644
index 63c88a34f6b7ec45193602838b2dc2568756f9ac..0000000000000000000000000000000000000000
--- a/Experiments/evaluate_naval.py
+++ /dev/null
@@ -1,81 +0,0 @@
-import os
-import numpy as np
-import torch
-import torch.backends.cudnn
-from torch.utils.data import DataLoader
-from torch.utils.tensorboard.writer import SummaryWriter
-
-from EIVArchitectures import Networks, initialize_weights
-from EIVData.naval_propulsion import load_data
-from EIVTrainingRoutines import train_and_store, loss_functions
-
-from train_noneiv_naval import p, init_std_y_list, seed_list, unscaled_reg, hidden_layers
-
-
-train_data, test_data = load_data()
-test_dataloader = DataLoader(test_data, batch_size=int(np.max((len(test_data), 800))))
-
-seed = seed_list[0]
-init_std_y = init_std_y_list[0]
-saved_file = os.path.join('saved_networks',
-            f'noneiv_naval'\
-                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
-                    f'_p_{p:.2f}_seed_{seed}.pkl')
-
-input_dim = train_data[0][0].numel()
-output_dim = train_data[0][1].numel()
-net = Networks.FNNBer(p=p, init_std_y=init_std_y,
-        h=[input_dim, *hidden_layers, output_dim])
-train_and_store.open_stored_training(saved_file=saved_file,
-        net=net)
-
-
-# # RMSE
-# x,y = next(iter(test_dataloader))
-# training_state = net.training
-# net.eval()
-# out = net(x)[0]
-# if len(y.shape) <=1:
-#     y = y.view((-1,1))
-# assert y.shape == out.shape
-# res = y-out
-# if training_state:
-#     net.train()
-# scale = train_data.dataset.std_labels
-# scaled_res = res * scale.view((1,-1))
-# scaled_res = scaled_res.detach().cpu().numpy().flatten()
-# rmse = np.sqrt(np.mean(scaled_res**2)) 
-# print(f'no Dropout RMSE {rmse:.3f}')
-
-# # RMSE with prediction
-# x,y = next(iter(test_dataloader))
-# training_state = net.training
-# net.train()
-# out = net.predict(x, number_of_draws=100)[0]
-# if len(y.shape) <=1:
-#     y = y.view((-1,1))
-# assert y.shape == out.shape
-# res = y-out
-# if training_state:
-#     net.train()
-# else:
-#     net.eval()
-# scale = train_data.dataset.std_labels
-# scaled_res = res * scale.view((1,-1))
-# scaled_res = scaled_res.detach().cpu().numpy().flatten()
-# rmse = np.sqrt(np.mean(scaled_res**2)) 
-# print(f'Dropout predictive RMSE {rmse:.3f}')
-
-
-# NLL
-x,y = next(iter(test_dataloader))
-training_state = net.training
-net.train()
-logdens = net.predictive_logdensity(x, y, number_of_draws=100,
-        decouple_dimensions=True,
-        scale_labels=train_data.dataset.std_labels.view((-1,))).mean()
-if training_state:
-    net.train()
-else:
-    net.eval()
-print(f'Dropout predictive {logdens:.3f}')
diff --git a/Experiments/evaluate_power.py b/Experiments/evaluate_power.py
deleted file mode 100644
index bc4d328af5cd624d5e2b56b66fd5dbbe4ba08d4a..0000000000000000000000000000000000000000
--- a/Experiments/evaluate_power.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import os
-import numpy as np
-import torch
-import torch.backends.cudnn
-from torch.utils.data import DataLoader
-from torch.utils.tensorboard.writer import SummaryWriter
-
-from EIVArchitectures import Networks, initialize_weights
-from EIVData.power_plant import load_data
-from EIVTrainingRoutines import train_and_store, loss_functions
-
-from train_noneiv_power import p, init_std_y_list, seed_list, unscaled_reg, hidden_layers
-
-
-train_data, test_data = load_data()
-test_dataloader = DataLoader(test_data, batch_size=int(np.max((len(test_data), 800))))
-
-seed = seed_list[0]
-init_std_y = init_std_y_list[0]
-saved_file = os.path.join('saved_networks',
-            f'noneiv_power'\
-                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
-                    f'_p_{p:.2f}_seed_{seed}.pkl')
-
-input_dim = train_data[0][0].numel()
-output_dim = train_data[0][1].numel()
-net = Networks.FNNBer(p=p, init_std_y=init_std_y,
-        h=[input_dim, *hidden_layers, output_dim])
-train_and_store.open_stored_training(saved_file=saved_file,
-        net=net)
-
-
-# RMSE
-x,y = next(iter(test_dataloader))
-out = net(x)[0]
-if len(y.shape) <=1:
-    y = y.view((-1,1))
-assert y.shape == out.shape
-res = y-out
-scale = train_data.dataset.std_labels
-scaled_res = res * scale.view((1,-1))
-scaled_res = scaled_res.detach().cpu().numpy().flatten()
-rmse = np.sqrt(np.mean(scaled_res**2)) 
-print(f'RMSE {rmse:.3f}')
-
-
-# NLL
-x,y = next(iter(test_dataloader))
-training_state = net.training
-net.train()
-logdens = net.predictive_logdensity(x, y, number_of_draws=100,
-        decouple_dimensions=True,
-        scale_labels=train_data.dataset.std_labels.view((-1,))).mean()
-if training_state:
-    net.train()
-else:
-    net.eval()
-print(f'Dropout predictive {logdens:.3f}')
diff --git a/Experiments/evaluate_protein.py b/Experiments/evaluate_protein.py
deleted file mode 100644
index de32d3cfbb99e949cda31e1451ec2e827adf4c66..0000000000000000000000000000000000000000
--- a/Experiments/evaluate_protein.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import os
-import numpy as np
-import torch
-import torch.backends.cudnn
-from torch.utils.data import DataLoader
-from torch.utils.tensorboard.writer import SummaryWriter
-
-from EIVArchitectures import Networks, initialize_weights
-from EIVData.protein_structure import load_data
-from EIVTrainingRoutines import train_and_store, loss_functions
-
-from train_noneiv_protein import p, init_std_y_list, seed_list, unscaled_reg, hidden_layers
-
-
-train_data, test_data = load_data()
-test_dataloader = DataLoader(test_data, batch_size=int(np.max((len(test_data), 800))))
-
-seed = seed_list[0]
-init_std_y = init_std_y_list[0]
-saved_file = os.path.join('saved_networks',
-            f'noneiv_protein'\
-                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
-                    f'_p_{p:.2f}_seed_{seed}.pkl')
-
-input_dim = train_data[0][0].numel()
-output_dim = train_data[0][1].numel()
-net = Networks.FNNBer(p=p, init_std_y=init_std_y,
-        h=[input_dim, *hidden_layers, output_dim])
-train_and_store.open_stored_training(saved_file=saved_file,
-        net=net)
-
-
-# RMSE
-x,y = next(iter(test_dataloader))
-out = net(x)[0]
-if len(y.shape) <=1:
-    y = y.view((-1,1))
-assert y.shape == out.shape
-res = y-out
-scale = train_data.dataset.std_labels
-scaled_res = res * scale.view((1,-1))
-scaled_res = scaled_res.detach().cpu().numpy().flatten()
-rmse = np.sqrt(np.mean(scaled_res**2)) 
-print(f'RMSE {rmse:.3f}')
-
-
-# NLL
-x,y = next(iter(test_dataloader))
-training_state = net.training
-net.train()
-logdens = net.predictive_logdensity(x, y, number_of_draws=100,
-        decouple_dimensions=True,
-        scale_labels=train_data.dataset.std_labels.view((-1,))).mean()
-if training_state:
-    net.train()
-else:
-    net.eval()
-print(f'Dropout predictive {logdens:.3f}')
diff --git a/Experiments/evaluate_tabular.py b/Experiments/evaluate_tabular.py
new file mode 100644
index 0000000000000000000000000000000000000000..d74c553789953779c12fee3a70f02de4b9d08328
--- /dev/null
+++ b/Experiments/evaluate_tabular.py
@@ -0,0 +1,191 @@
+import importlib
+import os
+import matplotlib
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from tqdm import tqdm
+
+from EIVArchitectures import Networks
+from EIVTrainingRoutines import train_and_store
+
+
+long_dataname = 'energy_efficiency'
+short_dataname = 'energy'
+
+load_data = importlib.import_module(f'EIVData.{long_dataname}').load_data
+train_noneiv = importlib.import_module(f'train_noneiv_{short_dataname}')
+train_eiv = importlib.import_module(f'train_eiv_{short_dataname}')
+
+train_data, test_data = load_data()
+test_dataloader = DataLoader(test_data, batch_size=int(np.max((len(test_data),
+    64))), shuffle=True)
+input_dim = train_data[0][0].numel()
+output_dim = train_data[0][1].numel()
+
+def collect_metrics(x,y, seed=0,
+    noneiv_number_of_draws=100, eiv_number_of_draws=[100,5],
+    decouple_dimensions=False, device=torch.device('cuda:1')):
+    """
+    :param x: A torch.tensor, taken as input
+    :param y: A torch.tensor, taken as output
+    :param seed: Integer. The seed used for loading, defaults to 0.
+    :param noneiv_number_of_draws: Number of draws for non-EiV model
+    for sampling from the posterior predictive. Defaults to 100.
+    :param noneiv_number_of_draws: Number of draws for EiV model
+    for sampling from the posterior predictive. Defaults to 100.
+    :param decouple_dimensions: Boolean. If True, the unsual convention
+    of Gal et al. is followed where, in the evaluation of the
+    log-posterior-predictive, each dimension is treated independently and then
+    averaged. If False (default), a multivariate distribution is used.
+    :returns: noneiv_rmse, noneiv_logdens, noneiv_bias,
+    eiv_rmse, eiv_logdens, eiv_bias
+    """
+    x,y = x.to(device), y.to(device)
+    init_std_y = train_noneiv.init_std_y_list[0]
+    unscaled_reg = train_noneiv.unscaled_reg
+    p = train_noneiv.p
+    hidden_layers = train_noneiv.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'_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)
+    train_and_store.open_stored_training(saved_file=saved_file,
+            net=net, device=device)
+
+
+    # RMSE
+    training_state = net.training
+    net.train()
+    out = net.predict(x, number_of_draws=noneiv_number_of_draws, 
+            take_average_of_prediction=True)[0]
+    if len(y.shape) <= 1:
+        y = y.view((-1,1))
+    assert y.shape == out.shape
+    res = y-out
+    scale = train_data.dataset.std_labels.to(device)
+    scaled_res = res * scale.view((1,-1))
+    scaled_res = scaled_res.detach().cpu().numpy().flatten()
+    noneiv_rmse = np.sqrt(np.mean(scaled_res**2))
+    noneiv_bias = np.mean(scaled_res)
+
+
+    # NLL
+    training_state = net.training
+    net.train()
+    noneiv_logdens = net.predictive_logdensity(x, y, number_of_draws=100,
+            decouple_dimensions=decouple_dimensions,
+            scale_labels=\
+                   train_data.dataset.std_labels.view((-1,)).to(device)\
+                   ).mean().detach().cpu().numpy()
+    if training_state:
+        net.train()
+    else:
+        net.eval()
+
+    # EiV
+    init_std_y = train_eiv.init_std_y_list[0]
+    unscaled_reg = train_eiv.unscaled_reg
+    p = train_eiv.p
+    hidden_layers = train_eiv.hidden_layers
+    fixed_std_x = train_eiv.fixed_std_x
+    saved_file = os.path.join('saved_networks',
+            f'eiv_{short_dataname}'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    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)
+    train_and_store.open_stored_training(saved_file=saved_file,
+            net=net)
+    # RMSE
+    training_state = net.training
+    noise_state = net.noise_is_on
+    net.train()
+    net.noise_on()
+    out = net.predict(x, number_of_draws=eiv_number_of_draws,
+            take_average_of_prediction=True)[0]
+    if len(y.shape) <=1:
+        y = y.view((-1,1))
+    assert y.shape == out.shape
+    res = y-out
+    scale = train_data.dataset.std_labels.to(device)
+    scaled_res = res * scale.view((1,-1))
+    scaled_res = scaled_res.detach().cpu().numpy().flatten()
+    eiv_rmse = np.sqrt(np.mean(scaled_res**2))
+    eiv_bias = np.mean(scaled_res)
+    if training_state:
+        net.train()
+    else:
+        net.eval()
+    if noise_state:
+        net.noise_on()
+    else:
+        net.noise_off()
+
+
+    # NLL
+    training_state = net.training
+    net.train()
+    eiv_logdens = net.predictive_logdensity(x, y,
+            number_of_draws=eiv_number_of_draws,
+            decouple_dimensions=decouple_dimensions,
+            scale_labels=\
+            train_data.dataset.std_labels.view((-1,)).to(device)\
+            ).mean().detach().cpu().numpy()
+    if training_state:
+        net.train()
+    else:
+        net.eval()
+    return noneiv_rmse, noneiv_logdens, noneiv_bias, \
+            eiv_rmse, eiv_logdens, eiv_bias
+
+noneiv_rmse_collection = []
+noneiv_logdens_collection = []
+noneiv_bias_collection = []
+eiv_rmse_collection = []
+eiv_logdens_collection = []
+eiv_bias_collection = []
+num_test_epochs = 10
+assert train_noneiv.seed_list == train_eiv.seed_list
+seed_list = train_noneiv.seed_list
+max_batch_number = 2
+for seed in tqdm(seed_list):
+    train_data, test_data = load_data(seed=seed)
+    test_dataloader = DataLoader(test_data,
+            batch_size=int(np.max((len(test_data),
+        800))), shuffle=True)
+    for i in tqdm(range(num_test_epochs)):
+        for j, (x,y) in enumerate(test_dataloader):
+            if j > max_batch_number:
+                break
+            noneiv_rmse, noneiv_logdens, noneiv_bias, \
+                    eiv_rmse, eiv_logdens, eiv_bias =\
+                    collect_metrics(x,y, seed=seed)
+            noneiv_rmse_collection.append(noneiv_rmse)
+            noneiv_logdens_collection.append(noneiv_logdens)
+            noneiv_bias_collection.append(noneiv_bias)
+            eiv_rmse_collection.append(eiv_rmse)
+            eiv_logdens_collection.append(eiv_logdens)
+            eiv_bias_collection.append(eiv_bias)
+
+
+print('Non-EiV')
+print(f'RMSE {np.mean(noneiv_rmse_collection):.5f}'\
+        f'({np.std(noneiv_rmse_collection)/np.sqrt(num_test_epochs*len(seed_list)):.5f})')
+print(f'LogDens {np.mean(noneiv_logdens_collection):.5f}'\
+        f'({np.std(noneiv_logdens_collection)/np.sqrt(num_test_epochs*len(seed_list)):.5f})')
+print(f'Bias {np.mean(noneiv_bias_collection):.5f}'\
+        f'({np.std(noneiv_bias_collection)/np.sqrt(num_test_epochs*len(seed_list)):.5f})')
+print('EiV')
+print(f'RMSE {np.mean(eiv_rmse_collection):.5f}'\
+        f'({np.std(eiv_rmse_collection)/np.sqrt(num_test_epochs*len(seed_list)):.5f})')
+print(f'LogDens {np.mean(eiv_logdens_collection):.5f}'\
+        f'({np.std(eiv_logdens_collection)/np.sqrt(num_test_epochs*len(seed_list)):.5f})')
+print(f'Bias {np.mean(eiv_bias_collection):.5f}'\
+        f'({np.std(eiv_bias_collection)/np.sqrt(num_test_epochs*len(seed_list)):.5f})')
diff --git a/Experiments/evaluate_wine.py b/Experiments/evaluate_wine.py
deleted file mode 100644
index 02be10de8141a299c02450b91deb87312aca4fb7..0000000000000000000000000000000000000000
--- a/Experiments/evaluate_wine.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import os
-import numpy as np
-import torch
-import torch.backends.cudnn
-from torch.utils.data import DataLoader
-from torch.utils.tensorboard.writer import SummaryWriter
-
-from EIVArchitectures import Networks, initialize_weights
-from EIVData.wine_quality import load_data
-from EIVTrainingRoutines import train_and_store, loss_functions
-
-from train_noneiv_wine import p, init_std_y_list, seed_list, unscaled_reg, hidden_layers
-
-
-train_data, test_data = load_data()
-test_dataloader = DataLoader(test_data, batch_size=int(np.max((len(test_data), 800))))
-
-seed = seed_list[0]
-init_std_y = init_std_y_list[0]
-saved_file = os.path.join('saved_networks',
-            f'noneiv_wine'\
-                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
-                    f'_p_{p:.2f}_seed_{seed}.pkl')
-
-input_dim = train_data[0][0].numel()
-output_dim = train_data[0][1].numel()
-net = Networks.FNNBer(p=p, init_std_y=init_std_y,
-        h=[input_dim, *hidden_layers, output_dim])
-train_and_store.open_stored_training(saved_file=saved_file,
-        net=net)
-
-
-# RMSE
-x,y = next(iter(test_dataloader))
-out = net(x)[0]
-if len(y.shape) <=1:
-    y = y.view((-1,1))
-assert y.shape == out.shape
-res = y-out
-scale = train_data.dataset.std_labels
-scaled_res = res * scale.view((1,-1))
-scaled_res = scaled_res.detach().cpu().numpy().flatten()
-rmse = np.sqrt(np.mean(scaled_res**2)) 
-print(f'RMSE {rmse:.3f}')
-
-
-# NLL
-x,y = next(iter(test_dataloader))
-training_state = net.training
-net.train()
-logdens = net.predictive_logdensity(x, y, number_of_draws=100,
-        decouple_dimensions=True,
-        scale_labels=train_data.dataset.std_labels.view((-1,))).mean()
-if training_state:
-    net.train()
-else:
-    net.eval()
-print(f'Dropout predictive {logdens:.3f}')
diff --git a/Experiments/evaluate_yacht.py b/Experiments/evaluate_yacht.py
deleted file mode 100644
index 842674dd5643d7032b32b487b0ac90c558121d72..0000000000000000000000000000000000000000
--- a/Experiments/evaluate_yacht.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import os
-import numpy as np
-import torch
-import torch.backends.cudnn
-from torch.utils.data import DataLoader
-from torch.utils.tensorboard.writer import SummaryWriter
-
-from EIVArchitectures import Networks, initialize_weights
-from EIVData.yacht_hydrodynamics import load_data
-from EIVTrainingRoutines import train_and_store, loss_functions
-
-from train_noneiv_yacht import p, init_std_y_list, seed_list, unscaled_reg, hidden_layers
-
-
-train_data, test_data = load_data()
-test_dataloader = DataLoader(test_data, batch_size=int(np.max((len(test_data), 800))))
-
-seed = seed_list[0]
-init_std_y = init_std_y_list[0]
-saved_file = os.path.join('saved_networks',
-            f'noneiv_yacht'\
-                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
-                    f'_p_{p:.2f}_seed_{seed}.pkl')
-
-input_dim = train_data[0][0].numel()
-output_dim = train_data[0][1].numel()
-net = Networks.FNNBer(p=p, init_std_y=init_std_y,
-        h=[input_dim, *hidden_layers, output_dim])
-train_and_store.open_stored_training(saved_file=saved_file,
-        net=net)
-
-
-# RMSE
-x,y = next(iter(test_dataloader))
-out = net(x)[0]
-if len(y.shape) <=1:
-    y = y.view((-1,1))
-assert y.shape == out.shape
-res = y-out
-scale = train_data.dataset.std_labels
-scaled_res = res * scale.view((1,-1))
-scaled_res = scaled_res.detach().cpu().numpy().flatten()
-rmse = np.sqrt(np.mean(scaled_res**2)) 
-print(f'RMSE {rmse:.3f}')
-
-
-# NLL
-x,y = next(iter(test_dataloader))
-training_state = net.training
-net.train()
-logdens = net.predictive_logdensity(x, y, number_of_draws=100,
-        decouple_dimensions=True,
-        scale_labels=train_data.dataset.std_labels.view((-1,))).mean()
-if training_state:
-    net.train()
-else:
-    net.eval()
-print(f'Dropout predictive {logdens:.3f}')
diff --git a/Experiments/train_eiv_california.py b/Experiments/train_eiv_california.py
new file mode 100644
index 0000000000000000000000000000000000000000..39c95fbe0cad3b7f839a8c230e9c54e5bfcde8f3
--- /dev/null
+++ b/Experiments/train_eiv_california.py
@@ -0,0 +1,155 @@
+"""
+Train EiV model on california housing dataset using different seeds
+"""
+import random
+import os
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from torch.utils.tensorboard.writer import SummaryWriter
+
+from EIVArchitectures import Networks, initialize_weights
+from EIVData.california_housing import load_data
+from EIVTrainingRoutines import train_and_store, loss_functions
+
+# hyperparameters
+lr = 1e-3
+batch_size = 200
+test_batch_size = 800
+number_of_epochs = 100
+unscaled_reg = 10
+report_point = 5
+p = 0.1
+lr_update = 20
+# pretraining = 300
+epoch_offset = 10
+init_std_y_list = [0.5]
+gamma = 0.5
+hidden_layers = [1024, 1024, 1024, 1024]
+device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
+fixed_std_x = 0.05
+
+# reproducability
+def set_seeds(seed):
+    torch.backends.cudnn.benchmark = False
+    np.random.seed(seed)
+    random.seed(seed) 
+    torch.manual_seed(seed)
+seed_list = range(10)
+
+# to store the RMSE
+rmse_chain = []
+
+class UpdatedTrainEpoch(train_and_store.TrainEpoch):
+    def pre_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch == 0:
+            self.lr = self.initial_lr
+            self.optimizer = torch.optim.Adam(net.parameters(), lr=self.lr)
+            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(
+            self.optimizer, lr_update, gamma)
+
+
+    def post_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch >= epoch_offset:
+            net.std_y_par.requires_grad = True
+        self.lr_scheduler.step() 
+
+    def extra_report(self, net, i):
+        """
+        Overwrites the corresponding method
+        and fed after initialization of this class
+        """
+        rmse = self.rmse(net).item()
+        rmse_chain.append(rmse)
+        writer.add_scalar('RMSE', rmse, self.total_count)
+        writer.add_scalar('train loss', self.last_train_loss, self.total_count)
+        writer.add_scalar('test loss', self.last_test_loss, self.total_count)
+        print(f'RMSE {rmse:.3f}')
+
+    def rmse(self, net):
+        """
+        Compute the root mean squared error for `net`
+        """
+        net_train_state = net.training
+        net_noise_state = net.noise_is_on
+        net.eval()
+        net.noise_off()
+        x, y = next(iter(self.test_dataloader))
+        if len(y.shape) <= 1:
+            y = y.view((-1,1))
+        out = net(x.to(device))[0].detach().cpu()
+        assert out.shape == y.shape
+        if net_train_state:
+            net.train()
+        if net_noise_state:
+            net.noise_on()
+        return torch.sqrt(torch.mean((out-y)**2))
+
+def train_on_data(init_std_y, seed):
+    """
+    Sets `seed`, loads data and trains an Bernoulli Modell, starting with
+    `init_std_y`.
+    """
+    # set seed
+    set_seeds(seed)
+    # load Datasets
+    train_data, test_data = load_data(seed=seed, splitting_part=0.8,
+            normalize=True)
+    # make dataloaders
+    train_dataloader = DataLoader(train_data, batch_size=batch_size, 
+            shuffle=True)
+    test_dataloader = DataLoader(test_data, batch_size=test_batch_size,
+            shuffle=True)
+    # create a net
+    input_dim = train_data[0][0].numel()
+    output_dim = train_data[0][1].numel()
+    net = Networks.FNNEIV(p=p,
+            init_std_y=init_std_y,
+            h=[input_dim, *hidden_layers, output_dim],
+            fixed_std_x=fixed_std_x)
+    net.apply(initialize_weights.glorot_init)
+    net = net.to(device)
+    net.std_y_par.requires_grad = False
+    std_x_map = lambda: net.get_std_x().detach().cpu().item()
+    std_y_map = lambda: net.get_std_y().detach().cpu().item()
+    # regularization
+    reg = unscaled_reg/len(train_data)
+    # create epoch_map
+    criterion = loss_functions.nll_eiv
+    epoch_map = UpdatedTrainEpoch(train_dataloader=train_dataloader,
+            test_dataloader=test_dataloader,
+            criterion=criterion, std_y_map=std_y_map, std_x_map=std_x_map,
+            lr=lr, reg=reg, report_point=report_point, device=device)
+    # run and save
+    save_file = os.path.join('saved_networks',
+            f'eiv_california'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    train_and_store.train_and_store(net=net, 
+            epoch_map=epoch_map,
+            number_of_epochs=number_of_epochs,
+            save_file=save_file)
+    
+
+if __name__ == '__main__':
+    for seed in seed_list:
+        # Tensorboard monitoring
+        writer = SummaryWriter(log_dir=f'/home/martin09/tmp/tensorboard/'\
+                f'run_eiv_california_lr_{lr:.4f}_seed'\
+                f'_{seed}_uregu_{unscaled_reg:.1f}_p_{p:.2f}'\
+                f'_fixed_std_x_{fixed_std_x:.3f}')
+        print(f'>>>>SEED: {seed}')
+        for init_std_y in init_std_y_list:
+            print(f'Using init_std_y={init_std_y:.3f}')
+            train_on_data(init_std_y, seed)
+
+
diff --git a/Experiments/train_eiv_concrete.py b/Experiments/train_eiv_concrete.py
new file mode 100644
index 0000000000000000000000000000000000000000..0cc42d5e05181c471b8281c2ff98eefc1a08aa41
--- /dev/null
+++ b/Experiments/train_eiv_concrete.py
@@ -0,0 +1,155 @@
+"""
+Train EiV model on concrete strength dataset using different seeds
+"""
+import random
+import os
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from torch.utils.tensorboard.writer import SummaryWriter
+
+from EIVArchitectures import Networks, initialize_weights
+from EIVData.concrete_strength import load_data
+from EIVTrainingRoutines import train_and_store, loss_functions
+
+# hyperparameters
+lr = 1e-3
+batch_size = 32
+test_batch_size = 800
+number_of_epochs = 100
+unscaled_reg = 10
+report_point = 5
+p = 0.2
+lr_update = 20
+# pretraining = 300
+epoch_offset = 10
+init_std_y_list = [0.5]
+gamma = 0.5
+hidden_layers = [1024, 1024, 1024, 1024]
+device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
+fixed_std_x = 0.05
+
+# reproducability
+def set_seeds(seed):
+    torch.backends.cudnn.benchmark = False
+    np.random.seed(seed)
+    random.seed(seed) 
+    torch.manual_seed(seed)
+seed_list = range(10)
+
+# to store the RMSE
+rmse_chain = []
+
+class UpdatedTrainEpoch(train_and_store.TrainEpoch):
+    def pre_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch == 0:
+            self.lr = self.initial_lr
+            self.optimizer = torch.optim.Adam(net.parameters(), lr=self.lr)
+            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(
+            self.optimizer, lr_update, gamma)
+
+
+    def post_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch >= epoch_offset:
+            net.std_y_par.requires_grad = True
+        self.lr_scheduler.step() 
+
+    def extra_report(self, net, i):
+        """
+        Overwrites the corresponding method
+        and fed after initialization of this class
+        """
+        rmse = self.rmse(net).item()
+        rmse_chain.append(rmse)
+        writer.add_scalar('RMSE', rmse, self.total_count)
+        writer.add_scalar('train loss', self.last_train_loss, self.total_count)
+        writer.add_scalar('test loss', self.last_test_loss, self.total_count)
+        print(f'RMSE {rmse:.3f}')
+
+    def rmse(self, net):
+        """
+        Compute the root mean squared error for `net`
+        """
+        net_train_state = net.training
+        net_noise_state = net.noise_is_on
+        net.eval()
+        net.noise_off()
+        x, y = next(iter(self.test_dataloader))
+        if len(y.shape) <= 1:
+            y = y.view((-1,1))
+        out = net(x.to(device))[0].detach().cpu()
+        assert out.shape == y.shape
+        if net_train_state:
+            net.train()
+        if net_noise_state:
+            net.noise_on()
+        return torch.sqrt(torch.mean((out-y)**2))
+
+def train_on_data(init_std_y, seed):
+    """
+    Sets `seed`, loads data and trains an Bernoulli Modell, starting with
+    `init_std_y`.
+    """
+    # set seed
+    set_seeds(seed)
+    # load Datasets
+    train_data, test_data = load_data(seed=seed, splitting_part=0.8,
+            normalize=True)
+    # make dataloaders
+    train_dataloader = DataLoader(train_data, batch_size=batch_size, 
+            shuffle=True)
+    test_dataloader = DataLoader(test_data, batch_size=test_batch_size,
+            shuffle=True)
+    # create a net
+    input_dim = train_data[0][0].numel()
+    output_dim = train_data[0][1].numel()
+    net = Networks.FNNEIV(p=p,
+            init_std_y=init_std_y,
+            h=[input_dim, *hidden_layers, output_dim],
+            fixed_std_x=fixed_std_x)
+    net.apply(initialize_weights.glorot_init)
+    net = net.to(device)
+    net.std_y_par.requires_grad = False
+    std_x_map = lambda: net.get_std_x().detach().cpu().item()
+    std_y_map = lambda: net.get_std_y().detach().cpu().item()
+    # regularization
+    reg = unscaled_reg/len(train_data)
+    # create epoch_map
+    criterion = loss_functions.nll_eiv
+    epoch_map = UpdatedTrainEpoch(train_dataloader=train_dataloader,
+            test_dataloader=test_dataloader,
+            criterion=criterion, std_y_map=std_y_map, std_x_map=std_x_map,
+            lr=lr, reg=reg, report_point=report_point, device=device)
+    # run and save
+    save_file = os.path.join('saved_networks',
+            f'eiv_concrete'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    train_and_store.train_and_store(net=net, 
+            epoch_map=epoch_map,
+            number_of_epochs=number_of_epochs,
+            save_file=save_file)
+    
+
+if __name__ == '__main__':
+    for seed in seed_list:
+        # Tensorboard monitoring
+        writer = SummaryWriter(log_dir=f'/home/martin09/tmp/tensorboard/'\
+                f'run_eiv_concrete_lr_{lr:.4f}_seed'\
+                f'_{seed}_uregu_{unscaled_reg:.1f}_p_{p:.2f}'\
+                f'_fixed_std_x_{fixed_std_x:.3f}')
+        print(f'>>>>SEED: {seed}')
+        for init_std_y in init_std_y_list:
+            print(f'Using init_std_y={init_std_y:.3f}')
+            train_on_data(init_std_y, seed)
+
+
diff --git a/Experiments/train_eiv_energy.py b/Experiments/train_eiv_energy.py
new file mode 100644
index 0000000000000000000000000000000000000000..40be4b283ed12849418e6c760618f79a70122bda
--- /dev/null
+++ b/Experiments/train_eiv_energy.py
@@ -0,0 +1,155 @@
+"""
+Train EiV model on the energy efficiency dataset using different seeds
+"""
+import random
+import os
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from torch.utils.tensorboard.writer import SummaryWriter
+
+from EIVArchitectures import Networks, initialize_weights
+from EIVData.energy_efficiency import load_data
+from EIVTrainingRoutines import train_and_store, loss_functions
+
+# hyperparameters
+lr = 1e-3
+batch_size = 32
+test_batch_size = 600
+number_of_epochs = 600
+unscaled_reg = 10
+report_point = 5
+p = 0.2
+lr_update = 100
+# pretraining = 300
+epoch_offset = 100
+init_std_y_list = [0.5]
+gamma = 0.5
+hidden_layers = [1024, 1024, 1024, 1024]
+device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
+fixed_std_x = 0.05
+
+# reproducability
+def set_seeds(seed):
+    torch.backends.cudnn.benchmark = False
+    np.random.seed(seed)
+    random.seed(seed) 
+    torch.manual_seed(seed)
+seed_list = range(10)
+
+# to store the RMSE
+rmse_chain = []
+
+class UpdatedTrainEpoch(train_and_store.TrainEpoch):
+    def pre_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch == 0:
+            self.lr = self.initial_lr
+            self.optimizer = torch.optim.Adam(net.parameters(), lr=self.lr)
+            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(
+            self.optimizer, lr_update, gamma)
+
+
+    def post_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch >= epoch_offset:
+            net.std_y_par.requires_grad = True
+        self.lr_scheduler.step() 
+
+    def extra_report(self, net, i):
+        """
+        Overwrites the corresponding method
+        and fed after initialization of this class
+        """
+        rmse = self.rmse(net).item()
+        rmse_chain.append(rmse)
+        writer.add_scalar('RMSE', rmse, self.total_count)
+        writer.add_scalar('train loss', self.last_train_loss, self.total_count)
+        writer.add_scalar('test loss', self.last_test_loss, self.total_count)
+        print(f'RMSE {rmse:.3f}')
+
+    def rmse(self, net):
+        """
+        Compute the root mean squared error for `net`
+        """
+        net_train_state = net.training
+        net_noise_state = net.noise_is_on
+        net.eval()
+        net.noise_off()
+        x, y = next(iter(self.test_dataloader))
+        if len(y.shape) <= 1:
+            y = y.view((-1,1))
+        out = net(x.to(device))[0].detach().cpu()
+        assert out.shape == y.shape
+        if net_train_state:
+            net.train()
+        if net_noise_state:
+            net.noise_on()
+        return torch.sqrt(torch.mean((out-y)**2))
+
+def train_on_data(init_std_y, seed):
+    """
+    Sets `seed`, loads data and trains an Bernoulli Modell, starting with
+    `init_std_y`.
+    """
+    # set seed
+    set_seeds(seed)
+    # load Datasets
+    train_data, test_data = load_data(seed=seed, splitting_part=0.8,
+            normalize=True)
+    # make dataloaders
+    train_dataloader = DataLoader(train_data, batch_size=batch_size, 
+            shuffle=True)
+    test_dataloader = DataLoader(test_data, batch_size=test_batch_size,
+            shuffle=True)
+    # create a net
+    input_dim = train_data[0][0].numel()
+    output_dim = train_data[0][1].numel()
+    net = Networks.FNNEIV(p=p,
+            init_std_y=init_std_y,
+            h=[input_dim, *hidden_layers, output_dim],
+            fixed_std_x=fixed_std_x)
+    net.apply(initialize_weights.glorot_init)
+    net = net.to(device)
+    net.std_y_par.requires_grad = False
+    std_x_map = lambda: net.get_std_x().detach().cpu().item()
+    std_y_map = lambda: net.get_std_y().detach().cpu().item()
+    # regularization
+    reg = unscaled_reg/len(train_data)
+    # create epoch_map
+    criterion = loss_functions.nll_eiv
+    epoch_map = UpdatedTrainEpoch(train_dataloader=train_dataloader,
+            test_dataloader=test_dataloader,
+            criterion=criterion, std_y_map=std_y_map, std_x_map=std_x_map,
+            lr=lr, reg=reg, report_point=report_point, device=device)
+    # run and save
+    save_file = os.path.join('saved_networks',
+            f'eiv_energy'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    train_and_store.train_and_store(net=net, 
+            epoch_map=epoch_map,
+            number_of_epochs=number_of_epochs,
+            save_file=save_file)
+    
+
+if __name__ == '__main__':
+    for seed in seed_list:
+        # Tensorboard monitoring
+        writer = SummaryWriter(log_dir=f'/home/martin09/tmp/tensorboard/'\
+                f'run_eiv_energy_lr_{lr:.4f}_seed'\
+                f'_{seed}_uregu_{unscaled_reg:.1f}_p_{p:.2f}'\
+                f'_fixed_std_x_{fixed_std_x:.3f}')
+        print(f'>>>>SEED: {seed}')
+        for init_std_y in init_std_y_list:
+            print(f'Using init_std_y={init_std_y:.3f}')
+            train_on_data(init_std_y, seed)
+
+
diff --git a/Experiments/train_eiv_kin8nm.py b/Experiments/train_eiv_kin8nm.py
new file mode 100644
index 0000000000000000000000000000000000000000..58d191c062f61eed6a4725bb7f364ed0303aeaaa
--- /dev/null
+++ b/Experiments/train_eiv_kin8nm.py
@@ -0,0 +1,155 @@
+"""
+Train EiV model on the kin8nm dataset using different seeds
+"""
+import random
+import os
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from torch.utils.tensorboard.writer import SummaryWriter
+
+from EIVArchitectures import Networks, initialize_weights
+from EIVData.kin8nm import load_data
+from EIVTrainingRoutines import train_and_store, loss_functions
+
+# hyperparameters
+lr = 1e-3
+batch_size = 32
+test_batch_size = 600
+number_of_epochs = 30
+unscaled_reg = 10
+report_point = 5
+p = 0.2
+lr_update = 20
+# pretraining = 300
+epoch_offset = 19
+init_std_y_list = [0.5]
+gamma = 0.5
+hidden_layers = [1024, 1024, 1024, 1024]
+device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
+fixed_std_x = 0.05
+
+# reproducability
+def set_seeds(seed):
+    torch.backends.cudnn.benchmark = False
+    np.random.seed(seed)
+    random.seed(seed) 
+    torch.manual_seed(seed)
+seed_list = range(10)
+
+# to store the RMSE
+rmse_chain = []
+
+class UpdatedTrainEpoch(train_and_store.TrainEpoch):
+    def pre_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch == 0:
+            self.lr = self.initial_lr
+            self.optimizer = torch.optim.Adam(net.parameters(), lr=self.lr)
+            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(
+            self.optimizer, lr_update, gamma)
+
+
+    def post_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch >= epoch_offset:
+            net.std_y_par.requires_grad = True
+        self.lr_scheduler.step() 
+
+    def extra_report(self, net, i):
+        """
+        Overwrites the corresponding method
+        and fed after initialization of this class
+        """
+        rmse = self.rmse(net).item()
+        rmse_chain.append(rmse)
+        writer.add_scalar('RMSE', rmse, self.total_count)
+        writer.add_scalar('train loss', self.last_train_loss, self.total_count)
+        writer.add_scalar('test loss', self.last_test_loss, self.total_count)
+        print(f'RMSE {rmse:.3f}')
+
+    def rmse(self, net):
+        """
+        Compute the root mean squared error for `net`
+        """
+        net_train_state = net.training
+        net_noise_state = net.noise_is_on
+        net.eval()
+        net.noise_off()
+        x, y = next(iter(self.test_dataloader))
+        if len(y.shape) <= 1:
+            y = y.view((-1,1))
+        out = net(x.to(device))[0].detach().cpu()
+        assert out.shape == y.shape
+        if net_train_state:
+            net.train()
+        if net_noise_state:
+            net.noise_on()
+        return torch.sqrt(torch.mean((out-y)**2))
+
+def train_on_data(init_std_y, seed):
+    """
+    Sets `seed`, loads data and trains an Bernoulli Modell, starting with
+    `init_std_y`.
+    """
+    # set seed
+    set_seeds(seed)
+    # load Datasets
+    train_data, test_data = load_data(seed=seed, splitting_part=0.8,
+            normalize=True)
+    # make dataloaders
+    train_dataloader = DataLoader(train_data, batch_size=batch_size, 
+            shuffle=True)
+    test_dataloader = DataLoader(test_data, batch_size=test_batch_size,
+            shuffle=True)
+    # create a net
+    input_dim = train_data[0][0].numel()
+    output_dim = train_data[0][1].numel()
+    net = Networks.FNNEIV(p=p,
+            init_std_y=init_std_y,
+            h=[input_dim, *hidden_layers, output_dim],
+            fixed_std_x=fixed_std_x)
+    net.apply(initialize_weights.glorot_init)
+    net = net.to(device)
+    net.std_y_par.requires_grad = False
+    std_x_map = lambda: net.get_std_x().detach().cpu().item()
+    std_y_map = lambda: net.get_std_y().detach().cpu().item()
+    # regularization
+    reg = unscaled_reg/len(train_data)
+    # create epoch_map
+    criterion = loss_functions.nll_eiv
+    epoch_map = UpdatedTrainEpoch(train_dataloader=train_dataloader,
+            test_dataloader=test_dataloader,
+            criterion=criterion, std_y_map=std_y_map, std_x_map=std_x_map,
+            lr=lr, reg=reg, report_point=report_point, device=device)
+    # run and save
+    save_file = os.path.join('saved_networks',
+            f'eiv_kin8nm'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    train_and_store.train_and_store(net=net, 
+            epoch_map=epoch_map,
+            number_of_epochs=number_of_epochs,
+            save_file=save_file)
+    
+
+if __name__ == '__main__':
+    for seed in seed_list:
+        # Tensorboard monitoring
+        writer = SummaryWriter(log_dir=f'/home/martin09/tmp/tensorboard/'\
+                f'run_eiv_kin8nm_lr_{lr:.4f}_seed'\
+                f'_{seed}_uregu_{unscaled_reg:.1f}_p_{p:.2f}'\
+                f'_fixed_std_x_{fixed_std_x:.3f}')
+        print(f'>>>>SEED: {seed}')
+        for init_std_y in init_std_y_list:
+            print(f'Using init_std_y={init_std_y:.3f}')
+            train_on_data(init_std_y, seed)
+
+
diff --git a/Experiments/train_eiv_msd.py b/Experiments/train_eiv_msd.py
new file mode 100644
index 0000000000000000000000000000000000000000..16e617f0a64fb8f9567ac2846a22cefff109b901
--- /dev/null
+++ b/Experiments/train_eiv_msd.py
@@ -0,0 +1,155 @@
+"""
+Train EiV model on the million song dataset using different seeds
+"""
+import random
+import os
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from torch.utils.tensorboard.writer import SummaryWriter
+
+from EIVArchitectures import Networks, initialize_weights
+from EIVData.million_song import load_data
+from EIVTrainingRoutines import train_and_store, loss_functions
+
+# hyperparameters
+lr = 1e-3
+batch_size = 100
+test_batch_size = 600
+number_of_epochs = 10
+unscaled_reg = 10
+report_point = 5
+p = 0.2
+lr_update = 4
+# pretraining = 300
+epoch_offset = 4
+init_std_y_list = [0.5]
+gamma = 0.5
+hidden_layers = [1024, 1024, 1024, 1024]
+device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
+fixed_std_x = 0.05
+
+# reproducability
+def set_seeds(seed):
+    torch.backends.cudnn.benchmark = False
+    np.random.seed(seed)
+    random.seed(seed) 
+    torch.manual_seed(seed)
+seed_list = range(10)
+
+# to store the RMSE
+rmse_chain = []
+
+class UpdatedTrainEpoch(train_and_store.TrainEpoch):
+    def pre_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch == 0:
+            self.lr = self.initial_lr
+            self.optimizer = torch.optim.Adam(net.parameters(), lr=self.lr)
+            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(
+            self.optimizer, lr_update, gamma)
+
+
+    def post_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch >= epoch_offset:
+            net.std_y_par.requires_grad = True
+        self.lr_scheduler.step() 
+
+    def extra_report(self, net, i):
+        """
+        Overwrites the corresponding method
+        and fed after initialization of this class
+        """
+        rmse = self.rmse(net).item()
+        rmse_chain.append(rmse)
+        writer.add_scalar('RMSE', rmse, self.total_count)
+        writer.add_scalar('train loss', self.last_train_loss, self.total_count)
+        writer.add_scalar('test loss', self.last_test_loss, self.total_count)
+        print(f'RMSE {rmse:.3f}')
+
+    def rmse(self, net):
+        """
+        Compute the root mean squared error for `net`
+        """
+        net_train_state = net.training
+        net_noise_state = net.noise_is_on
+        net.eval()
+        net.noise_off()
+        x, y = next(iter(self.test_dataloader))
+        if len(y.shape) <= 1:
+            y = y.view((-1,1))
+        out = net(x.to(device))[0].detach().cpu()
+        assert out.shape == y.shape
+        if net_train_state:
+            net.train()
+        if net_noise_state:
+            net.noise_on()
+        return torch.sqrt(torch.mean((out-y)**2))
+
+def train_on_data(init_std_y, seed):
+    """
+    Sets `seed`, loads data and trains an Bernoulli Modell, starting with
+    `init_std_y`.
+    """
+    # set seed
+    set_seeds(seed)
+    # load Datasets
+    train_data, test_data = load_data(seed=seed, splitting_part=0.8,
+            normalize=True)
+    # make dataloaders
+    train_dataloader = DataLoader(train_data, batch_size=batch_size, 
+            shuffle=True)
+    test_dataloader = DataLoader(test_data, batch_size=test_batch_size,
+            shuffle=True)
+    # create a net
+    input_dim = train_data[0][0].numel()
+    output_dim = train_data[0][1].numel()
+    net = Networks.FNNEIV(p=p,
+            init_std_y=init_std_y,
+            h=[input_dim, *hidden_layers, output_dim],
+            fixed_std_x=fixed_std_x)
+    net.apply(initialize_weights.glorot_init)
+    net = net.to(device)
+    net.std_y_par.requires_grad = False
+    std_x_map = lambda: net.get_std_x().detach().cpu().item()
+    std_y_map = lambda: net.get_std_y().detach().cpu().item()
+    # regularization
+    reg = unscaled_reg/len(train_data)
+    # create epoch_map
+    criterion = loss_functions.nll_eiv
+    epoch_map = UpdatedTrainEpoch(train_dataloader=train_dataloader,
+            test_dataloader=test_dataloader,
+            criterion=criterion, std_y_map=std_y_map, std_x_map=std_x_map,
+            lr=lr, reg=reg, report_point=report_point, device=device)
+    # run and save
+    save_file = os.path.join('saved_networks',
+            f'eiv_msd'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    train_and_store.train_and_store(net=net, 
+            epoch_map=epoch_map,
+            number_of_epochs=number_of_epochs,
+            save_file=save_file)
+    
+
+if __name__ == '__main__':
+    for seed in seed_list:
+        # Tensorboard monitoring
+        writer = SummaryWriter(log_dir=f'/home/martin09/tmp/tensorboard/'\
+                f'run_eiv_msd_lr_{lr:.4f}_seed'\
+                f'_{seed}_uregu_{unscaled_reg:.1f}_p_{p:.2f}'\
+                f'_fixed_std_x_{fixed_std_x:.3f}')
+        print(f'>>>>SEED: {seed}')
+        for init_std_y in init_std_y_list:
+            print(f'Using init_std_y={init_std_y:.3f}')
+            train_on_data(init_std_y, seed)
+
+
diff --git a/Experiments/train_eiv_naval.py b/Experiments/train_eiv_naval.py
new file mode 100644
index 0000000000000000000000000000000000000000..8db266c735d022fe4a2170a7de4d0a5ba3610d74
--- /dev/null
+++ b/Experiments/train_eiv_naval.py
@@ -0,0 +1,155 @@
+"""
+Train EiV model on the naval propulsion dataset using different seeds
+"""
+import random
+import os
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from torch.utils.tensorboard.writer import SummaryWriter
+
+from EIVArchitectures import Networks, initialize_weights
+from EIVData.naval_propulsion import load_data
+from EIVTrainingRoutines import train_and_store, loss_functions
+
+# hyperparameters
+lr = 1e-3
+batch_size = 32
+test_batch_size = 600
+number_of_epochs = 30
+unscaled_reg = 10
+report_point = 5
+p = 0.2
+lr_update = 20
+# pretraining = 300
+epoch_offset = 20
+init_std_y_list = [0.5]
+gamma = 0.5
+hidden_layers = [1024, 1024, 1024, 1024]
+device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
+fixed_std_x = 0.05
+
+# reproducability
+def set_seeds(seed):
+    torch.backends.cudnn.benchmark = False
+    np.random.seed(seed)
+    random.seed(seed) 
+    torch.manual_seed(seed)
+seed_list = range(10)
+
+# to store the RMSE
+rmse_chain = []
+
+class UpdatedTrainEpoch(train_and_store.TrainEpoch):
+    def pre_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch == 0:
+            self.lr = self.initial_lr
+            self.optimizer = torch.optim.Adam(net.parameters(), lr=self.lr)
+            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(
+            self.optimizer, lr_update, gamma)
+
+
+    def post_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch >= epoch_offset:
+            net.std_y_par.requires_grad = True
+        self.lr_scheduler.step() 
+
+    def extra_report(self, net, i):
+        """
+        Overwrites the corresponding method
+        and fed after initialization of this class
+        """
+        rmse = self.rmse(net).item()
+        rmse_chain.append(rmse)
+        writer.add_scalar('RMSE', rmse, self.total_count)
+        writer.add_scalar('train loss', self.last_train_loss, self.total_count)
+        writer.add_scalar('test loss', self.last_test_loss, self.total_count)
+        print(f'RMSE {rmse:.3f}')
+
+    def rmse(self, net):
+        """
+        Compute the root mean squared error for `net`
+        """
+        net_train_state = net.training
+        net_noise_state = net.noise_is_on
+        net.eval()
+        net.noise_off()
+        x, y = next(iter(self.test_dataloader))
+        if len(y.shape) <= 1:
+            y = y.view((-1,1))
+        out = net(x.to(device))[0].detach().cpu()
+        assert out.shape == y.shape
+        if net_train_state:
+            net.train()
+        if net_noise_state:
+            net.noise_on()
+        return torch.sqrt(torch.mean((out-y)**2))
+
+def train_on_data(init_std_y, seed):
+    """
+    Sets `seed`, loads data and trains an Bernoulli Modell, starting with
+    `init_std_y`.
+    """
+    # set seed
+    set_seeds(seed)
+    # load Datasets
+    train_data, test_data = load_data(seed=seed, splitting_part=0.8,
+            normalize=True)
+    # make dataloaders
+    train_dataloader = DataLoader(train_data, batch_size=batch_size, 
+            shuffle=True)
+    test_dataloader = DataLoader(test_data, batch_size=test_batch_size,
+            shuffle=True)
+    # create a net
+    input_dim = train_data[0][0].numel()
+    output_dim = train_data[0][1].numel()
+    net = Networks.FNNEIV(p=p,
+            init_std_y=init_std_y,
+            h=[input_dim, *hidden_layers, output_dim],
+            fixed_std_x=fixed_std_x)
+    net.apply(initialize_weights.glorot_init)
+    net = net.to(device)
+    net.std_y_par.requires_grad = False
+    std_x_map = lambda: net.get_std_x().detach().cpu().item()
+    std_y_map = lambda: net.get_std_y().detach().cpu().item()
+    # regularization
+    reg = unscaled_reg/len(train_data)
+    # create epoch_map
+    criterion = loss_functions.nll_eiv
+    epoch_map = UpdatedTrainEpoch(train_dataloader=train_dataloader,
+            test_dataloader=test_dataloader,
+            criterion=criterion, std_y_map=std_y_map, std_x_map=std_x_map,
+            lr=lr, reg=reg, report_point=report_point, device=device)
+    # run and save
+    save_file = os.path.join('saved_networks',
+            f'eiv_naval'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    train_and_store.train_and_store(net=net, 
+            epoch_map=epoch_map,
+            number_of_epochs=number_of_epochs,
+            save_file=save_file)
+    
+
+if __name__ == '__main__':
+    for seed in seed_list:
+        # Tensorboard monitoring
+        writer = SummaryWriter(log_dir=f'/home/martin09/tmp/tensorboard/'\
+                f'run_eiv_naval_lr_{lr:.4f}_seed'\
+                f'_{seed}_uregu_{unscaled_reg:.1f}_p_{p:.2f}'\
+                f'_fixed_std_x_{fixed_std_x:.3f}')
+        print(f'>>>>SEED: {seed}')
+        for init_std_y in init_std_y_list:
+            print(f'Using init_std_y={init_std_y:.3f}')
+            train_on_data(init_std_y, seed)
+
+
diff --git a/Experiments/train_eiv_power.py b/Experiments/train_eiv_power.py
new file mode 100644
index 0000000000000000000000000000000000000000..2ddf75af75db7607a21cef6e212c9601d391b421
--- /dev/null
+++ b/Experiments/train_eiv_power.py
@@ -0,0 +1,153 @@
+"""
+Train EiV model on power plant dataset using different seeds
+"""
+import random
+import os
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from torch.utils.tensorboard.writer import SummaryWriter
+
+from EIVArchitectures import Networks, initialize_weights
+from EIVData.power_plant import load_data
+from EIVTrainingRoutines import train_and_store, loss_functions
+
+# hyperparameters
+lr = 1e-3
+batch_size = 64
+test_batch_size = 600
+number_of_epochs = 35
+unscaled_reg = 10
+report_point = 5
+p = 0.2
+lr_update = 10
+# pretraining = 300
+epoch_offset = 15
+init_std_y_list = [0.5]
+gamma = 0.5
+hidden_layers = [1024, 1024, 1024, 1024]
+device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
+fixed_std_x = 0.05
+
+# reproducability
+def set_seeds(seed):
+    torch.backends.cudnn.benchmark = False
+    np.random.seed(seed)
+    random.seed(seed) 
+    torch.manual_seed(seed)
+seed_list = range(10)
+
+# to store the RMSE
+rmse_chain = []
+
+class UpdatedTrainEpoch(train_and_store.TrainEpoch):
+    def pre_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch == 0:
+            self.lr = self.initial_lr
+            self.optimizer = torch.optim.Adam(net.parameters(), lr=self.lr)
+            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(
+            self.optimizer, lr_update, gamma)
+
+
+    def post_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch >= epoch_offset:
+            net.std_y_par.requires_grad = True
+        self.lr_scheduler.step() 
+
+    def extra_report(self, net, i):
+        """
+        Overwrites the corresponding method
+        and fed after initialization of this class
+        """
+        rmse = self.rmse(net).item()
+        rmse_chain.append(rmse)
+        writer.add_scalar('RMSE', rmse, self.total_count)
+        writer.add_scalar('train loss', self.last_train_loss, self.total_count)
+        writer.add_scalar('test loss', self.last_test_loss, self.total_count)
+        print(f'RMSE {rmse:.3f}')
+
+    def rmse(self, net):
+        """
+        Compute the root mean squared error for `net`
+        """
+        net_train_state = net.training
+        net_noise_state = net.noise_is_on
+        net.eval()
+        net.noise_off()
+        x, y = next(iter(self.test_dataloader))
+        if len(y.shape) <= 1:
+            y = y.view((-1,1))
+        out = net(x.to(device))[0].detach().cpu()
+        assert out.shape == y.shape
+        if net_train_state:
+            net.train()
+        if net_noise_state:
+            net.noise_on()
+        return torch.sqrt(torch.mean((out-y)**2))
+
+def train_on_data(init_std_y, seed):
+    """
+    Sets `seed`, loads data and trains an Bernoulli Modell, starting with
+    `init_std_y`.
+    """
+    # set seed
+    set_seeds(seed)
+    # load Datasets
+    train_data, test_data = load_data(seed=seed, splitting_part=0.8,
+            normalize=True)
+    # make dataloaders
+    train_dataloader = DataLoader(train_data, batch_size=batch_size, 
+            shuffle=True)
+    test_dataloader = DataLoader(test_data, batch_size=test_batch_size,
+            shuffle=True)
+    # create a net
+    input_dim = train_data[0][0].numel()
+    output_dim = train_data[0][1].numel()
+    net = Networks.FNNEIV(p=p,
+            init_std_y=init_std_y,
+            h=[input_dim, *hidden_layers, output_dim],
+            fixed_std_x=fixed_std_x)
+    net.apply(initialize_weights.glorot_init)
+    net = net.to(device)
+    net.std_y_par.requires_grad = False
+    std_x_map = lambda: net.get_std_x().detach().cpu().item()
+    std_y_map = lambda: net.get_std_y().detach().cpu().item()
+    # regularization
+    reg = unscaled_reg/len(train_data)
+    # create epoch_map
+    criterion = loss_functions.nll_eiv
+    epoch_map = UpdatedTrainEpoch(train_dataloader=train_dataloader,
+            test_dataloader=test_dataloader,
+            criterion=criterion, std_y_map=std_y_map, std_x_map=std_x_map,
+            lr=lr, reg=reg, report_point=report_point, device=device)
+    # run and save
+    save_file = os.path.join('saved_networks',
+            f'eiv_power'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    train_and_store.train_and_store(net=net, 
+            epoch_map=epoch_map,
+            number_of_epochs=number_of_epochs,
+            save_file=save_file)
+    
+
+if __name__ == '__main__':
+    for seed in seed_list:
+        # Tensorboard monitoring
+        writer = SummaryWriter(log_dir=f'/home/martin09/tmp/tensorboard/'\
+                f'run_eiv_power_lr_{lr:.4f}_seed'\
+                f'_{seed}_uregu_{unscaled_reg:.1f}_p_{p:.2f}'\
+                f'_fixed_std_x_{fixed_std_x:.3f}')
+        print(f'>>>>SEED: {seed}')
+        for init_std_y in init_std_y_list:
+            print(f'Using init_std_y={init_std_y:.3f}')
+            train_on_data(init_std_y, seed)
diff --git a/Experiments/train_eiv_protein.py b/Experiments/train_eiv_protein.py
new file mode 100644
index 0000000000000000000000000000000000000000..3801d5e07572a7f37b23b457a3adc8b83f3f8418
--- /dev/null
+++ b/Experiments/train_eiv_protein.py
@@ -0,0 +1,155 @@
+"""
+Train EiV model on protein structure dataset using different seeds
+"""
+import random
+import os
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from torch.utils.tensorboard.writer import SummaryWriter
+
+from EIVArchitectures import Networks, initialize_weights
+from EIVData.protein_structure import load_data
+from EIVTrainingRoutines import train_and_store, loss_functions
+
+# hyperparameters
+lr = 1e-3
+batch_size = 100
+test_batch_size = 600
+number_of_epochs = 30
+unscaled_reg = 10
+report_point = 5
+p = 0.2
+lr_update = 10
+# pretraining = 300
+epoch_offset = 10
+init_std_y_list = [0.5]
+gamma = 0.5
+hidden_layers = [1024, 1024, 1024, 1024]
+device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
+fixed_std_x = 0.05
+
+# reproducability
+def set_seeds(seed):
+    torch.backends.cudnn.benchmark = False
+    np.random.seed(seed)
+    random.seed(seed) 
+    torch.manual_seed(seed)
+seed_list = range(10)
+
+# to store the RMSE
+rmse_chain = []
+
+class UpdatedTrainEpoch(train_and_store.TrainEpoch):
+    def pre_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch == 0:
+            self.lr = self.initial_lr
+            self.optimizer = torch.optim.Adam(net.parameters(), lr=self.lr)
+            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(
+            self.optimizer, lr_update, gamma)
+
+
+    def post_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch >= epoch_offset:
+            net.std_y_par.requires_grad = True
+        self.lr_scheduler.step() 
+
+    def extra_report(self, net, i):
+        """
+        Overwrites the corresponding method
+        and fed after initialization of this class
+        """
+        rmse = self.rmse(net).item()
+        rmse_chain.append(rmse)
+        writer.add_scalar('RMSE', rmse, self.total_count)
+        writer.add_scalar('train loss', self.last_train_loss, self.total_count)
+        writer.add_scalar('test loss', self.last_test_loss, self.total_count)
+        print(f'RMSE {rmse:.3f}')
+
+    def rmse(self, net):
+        """
+        Compute the root mean squared error for `net`
+        """
+        net_train_state = net.training
+        net_noise_state = net.noise_is_on
+        net.eval()
+        net.noise_off()
+        x, y = next(iter(self.test_dataloader))
+        if len(y.shape) <= 1:
+            y = y.view((-1,1))
+        out = net(x.to(device))[0].detach().cpu()
+        assert out.shape == y.shape
+        if net_train_state:
+            net.train()
+        if net_noise_state:
+            net.noise_on()
+        return torch.sqrt(torch.mean((out-y)**2))
+
+def train_on_data(init_std_y, seed):
+    """
+    Sets `seed`, loads data and trains an Bernoulli Modell, starting with
+    `init_std_y`.
+    """
+    # set seed
+    set_seeds(seed)
+    # load Datasets
+    train_data, test_data = load_data(seed=seed, splitting_part=0.8,
+            normalize=True)
+    # make dataloaders
+    train_dataloader = DataLoader(train_data, batch_size=batch_size, 
+            shuffle=True)
+    test_dataloader = DataLoader(test_data, batch_size=test_batch_size,
+            shuffle=True)
+    # create a net
+    input_dim = train_data[0][0].numel()
+    output_dim = train_data[0][1].numel()
+    net = Networks.FNNEIV(p=p,
+            init_std_y=init_std_y,
+            h=[input_dim, *hidden_layers, output_dim],
+            fixed_std_x=fixed_std_x)
+    net.apply(initialize_weights.glorot_init)
+    net = net.to(device)
+    net.std_y_par.requires_grad = False
+    std_x_map = lambda: net.get_std_x().detach().cpu().item()
+    std_y_map = lambda: net.get_std_y().detach().cpu().item()
+    # regularization
+    reg = unscaled_reg/len(train_data)
+    # create epoch_map
+    criterion = loss_functions.nll_eiv
+    epoch_map = UpdatedTrainEpoch(train_dataloader=train_dataloader,
+            test_dataloader=test_dataloader,
+            criterion=criterion, std_y_map=std_y_map, std_x_map=std_x_map,
+            lr=lr, reg=reg, report_point=report_point, device=device)
+    # run and save
+    save_file = os.path.join('saved_networks',
+            f'eiv_protein'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    train_and_store.train_and_store(net=net, 
+            epoch_map=epoch_map,
+            number_of_epochs=number_of_epochs,
+            save_file=save_file)
+    
+
+if __name__ == '__main__':
+    for seed in seed_list:
+        # Tensorboard monitoring
+        writer = SummaryWriter(log_dir=f'/home/martin09/tmp/tensorboard/'\
+                f'run_eiv_protein_lr_{lr:.4f}_seed'\
+                f'_{seed}_uregu_{unscaled_reg:.1f}_p_{p:.2f}'\
+                f'_fixed_std_x_{fixed_std_x:.3f}')
+        print(f'>>>>SEED: {seed}')
+        for init_std_y in init_std_y_list:
+            print(f'Using init_std_y={init_std_y:.3f}')
+            train_on_data(init_std_y, seed)
+
+
diff --git a/Experiments/train_eiv_wine.py b/Experiments/train_eiv_wine.py
new file mode 100644
index 0000000000000000000000000000000000000000..ca561adc2e4d3f9a37060d90b370d21a06ed01d0
--- /dev/null
+++ b/Experiments/train_eiv_wine.py
@@ -0,0 +1,155 @@
+"""
+Train EiV model on wine quality dataset using different seeds
+"""
+import random
+import os
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from torch.utils.tensorboard.writer import SummaryWriter
+
+from EIVArchitectures import Networks, initialize_weights
+from EIVData.wine_quality import load_data
+from EIVTrainingRoutines import train_and_store, loss_functions
+
+# hyperparameters
+lr = 1e-3
+batch_size = 32
+test_batch_size = 800
+number_of_epochs = 100
+unscaled_reg = 10
+report_point = 5
+p = 0.2
+lr_update = 30
+# pretraining = 300
+epoch_offset = 50
+init_std_y_list = [0.5]
+gamma = 0.5
+hidden_layers = [1024, 1024, 1024, 1024]
+device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
+fixed_std_x = 0.05
+
+# reproducability
+def set_seeds(seed):
+    torch.backends.cudnn.benchmark = False
+    np.random.seed(seed)
+    random.seed(seed) 
+    torch.manual_seed(seed)
+seed_list = range(10)
+
+# to store the RMSE
+rmse_chain = []
+
+class UpdatedTrainEpoch(train_and_store.TrainEpoch):
+    def pre_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch == 0:
+            self.lr = self.initial_lr
+            self.optimizer = torch.optim.Adam(net.parameters(), lr=self.lr)
+            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(
+            self.optimizer, lr_update, gamma)
+
+
+    def post_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch >= epoch_offset:
+            net.std_y_par.requires_grad = True
+        self.lr_scheduler.step() 
+
+    def extra_report(self, net, i):
+        """
+        Overwrites the corresponding method
+        and fed after initialization of this class
+        """
+        rmse = self.rmse(net).item()
+        rmse_chain.append(rmse)
+        writer.add_scalar('RMSE', rmse, self.total_count)
+        writer.add_scalar('train loss', self.last_train_loss, self.total_count)
+        writer.add_scalar('test loss', self.last_test_loss, self.total_count)
+        print(f'RMSE {rmse:.3f}')
+
+    def rmse(self, net):
+        """
+        Compute the root mean squared error for `net`
+        """
+        net_train_state = net.training
+        net_noise_state = net.noise_is_on
+        net.eval()
+        net.noise_off()
+        x, y = next(iter(self.test_dataloader))
+        if len(y.shape) <= 1:
+            y = y.view((-1,1))
+        out = net(x.to(device))[0].detach().cpu()
+        assert out.shape == y.shape
+        if net_train_state:
+            net.train()
+        if net_noise_state:
+            net.noise_on()
+        return torch.sqrt(torch.mean((out-y)**2))
+
+def train_on_data(init_std_y, seed):
+    """
+    Sets `seed`, loads data and trains an Bernoulli Modell, starting with
+    `init_std_y`.
+    """
+    # set seed
+    set_seeds(seed)
+    # load Datasets
+    train_data, test_data = load_data(seed=seed, splitting_part=0.8,
+            normalize=True)
+    # make dataloaders
+    train_dataloader = DataLoader(train_data, batch_size=batch_size, 
+            shuffle=True)
+    test_dataloader = DataLoader(test_data, batch_size=test_batch_size,
+            shuffle=True)
+    # create a net
+    input_dim = train_data[0][0].numel()
+    output_dim = train_data[0][1].numel()
+    net = Networks.FNNEIV(p=p,
+            init_std_y=init_std_y,
+            h=[input_dim, *hidden_layers, output_dim],
+            fixed_std_x=fixed_std_x)
+    net.apply(initialize_weights.glorot_init)
+    net = net.to(device)
+    net.std_y_par.requires_grad = False
+    std_x_map = lambda: net.get_std_x().detach().cpu().item()
+    std_y_map = lambda: net.get_std_y().detach().cpu().item()
+    # regularization
+    reg = unscaled_reg/len(train_data)
+    # create epoch_map
+    criterion = loss_functions.nll_eiv
+    epoch_map = UpdatedTrainEpoch(train_dataloader=train_dataloader,
+            test_dataloader=test_dataloader,
+            criterion=criterion, std_y_map=std_y_map, std_x_map=std_x_map,
+            lr=lr, reg=reg, report_point=report_point, device=device)
+    # run and save
+    save_file = os.path.join('saved_networks',
+            f'eiv_wine'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    train_and_store.train_and_store(net=net, 
+            epoch_map=epoch_map,
+            number_of_epochs=number_of_epochs,
+            save_file=save_file)
+    
+
+if __name__ == '__main__':
+    for seed in seed_list:
+        # Tensorboard monitoring
+        writer = SummaryWriter(log_dir=f'/home/martin09/tmp/tensorboard/'\
+                f'run_eiv_wine_lr_{lr:.4f}_seed'\
+                f'_{seed}_uregu_{unscaled_reg:.1f}_p_{p:.2f}'\
+                f'_fixed_std_x_{fixed_std_x:.3f}')
+        print(f'>>>>SEED: {seed}')
+        for init_std_y in init_std_y_list:
+            print(f'Using init_std_y={init_std_y:.3f}')
+            train_on_data(init_std_y, seed)
+
+
diff --git a/Experiments/train_eiv_yacht.py b/Experiments/train_eiv_yacht.py
new file mode 100644
index 0000000000000000000000000000000000000000..2c9ee88853d77b99fefc3f2507ad827597895d61
--- /dev/null
+++ b/Experiments/train_eiv_yacht.py
@@ -0,0 +1,153 @@
+"""
+Train EiV model on the yacht hydrodynamics dataset using different seeds
+"""
+import random
+import os
+
+import numpy as np
+import torch
+import torch.backends.cudnn
+from torch.utils.data import DataLoader
+from torch.utils.tensorboard.writer import SummaryWriter
+
+from EIVArchitectures import Networks, initialize_weights
+from EIVData.yacht_hydrodynamics import load_data
+from EIVTrainingRoutines import train_and_store, loss_functions
+
+# hyperparameters
+lr = 1e-3
+batch_size = 32
+test_batch_size = 600
+number_of_epochs = 1200
+unscaled_reg = 10
+report_point = 5
+p = 0.2
+lr_update = 200
+# pretraining = 300
+epoch_offset = 250
+init_std_y_list = [0.5]
+gamma = 0.5
+hidden_layers = [1024, 1024, 1024, 1024]
+device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
+fixed_std_x = 0.05
+
+# reproducability
+def set_seeds(seed):
+    torch.backends.cudnn.benchmark = False
+    np.random.seed(seed)
+    random.seed(seed) 
+    torch.manual_seed(seed)
+seed_list = range(10)
+
+# to store the RMSE
+rmse_chain = []
+
+class UpdatedTrainEpoch(train_and_store.TrainEpoch):
+    def pre_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch == 0:
+            self.lr = self.initial_lr
+            self.optimizer = torch.optim.Adam(net.parameters(), lr=self.lr)
+            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(
+            self.optimizer, lr_update, gamma)
+
+
+    def post_epoch_update(self, net, epoch):
+        """
+        Overwrites the corresponding method
+        """
+        if epoch >= epoch_offset:
+            net.std_y_par.requires_grad = True
+        self.lr_scheduler.step() 
+
+    def extra_report(self, net, i):
+        """
+        Overwrites the corresponding method
+        and fed after initialization of this class
+        """
+        rmse = self.rmse(net).item()
+        rmse_chain.append(rmse)
+        writer.add_scalar('RMSE', rmse, self.total_count)
+        writer.add_scalar('train loss', self.last_train_loss, self.total_count)
+        writer.add_scalar('test loss', self.last_test_loss, self.total_count)
+        print(f'RMSE {rmse:.3f}')
+
+    def rmse(self, net):
+        """
+        Compute the root mean squared error for `net`
+        """
+        net_train_state = net.training
+        net_noise_state = net.noise_is_on
+        net.eval()
+        net.noise_off()
+        x, y = next(iter(self.test_dataloader))
+        if len(y.shape) <= 1:
+            y = y.view((-1,1))
+        out = net(x.to(device))[0].detach().cpu()
+        assert out.shape == y.shape
+        if net_train_state:
+            net.train()
+        if net_noise_state:
+            net.noise_on()
+        return torch.sqrt(torch.mean((out-y)**2))
+
+def train_on_data(init_std_y, seed):
+    """
+    Sets `seed`, loads data and trains an Bernoulli Modell, starting with
+    `init_std_y`.
+    """
+    # set seed
+    set_seeds(seed)
+    # load Datasets
+    train_data, test_data = load_data(seed=seed, splitting_part=0.8,
+            normalize=True)
+    # make dataloaders
+    train_dataloader = DataLoader(train_data, batch_size=batch_size, 
+            shuffle=True)
+    test_dataloader = DataLoader(test_data, batch_size=test_batch_size,
+            shuffle=True)
+    # create a net
+    input_dim = train_data[0][0].numel()
+    output_dim = train_data[0][1].numel()
+    net = Networks.FNNEIV(p=p,
+            init_std_y=init_std_y,
+            h=[input_dim, *hidden_layers, output_dim],
+            fixed_std_x=fixed_std_x)
+    net.apply(initialize_weights.glorot_init)
+    net = net.to(device)
+    net.std_y_par.requires_grad = False
+    std_x_map = lambda: net.get_std_x().detach().cpu().item()
+    std_y_map = lambda: net.get_std_y().detach().cpu().item()
+    # regularization
+    reg = unscaled_reg/len(train_data)
+    # create epoch_map
+    criterion = loss_functions.nll_eiv
+    epoch_map = UpdatedTrainEpoch(train_dataloader=train_dataloader,
+            test_dataloader=test_dataloader,
+            criterion=criterion, std_y_map=std_y_map, std_x_map=std_x_map,
+            lr=lr, reg=reg, report_point=report_point, device=device)
+    # run and save
+    save_file = os.path.join('saved_networks',
+            f'eiv_yacht'\
+                    f'_init_std_y_{init_std_y:.3f}_ureg_{unscaled_reg:.1f}'\
+                    f'_p_{p:.2f}_fixed_std_x_{fixed_std_x:.3f}'\
+                    f'_seed_{seed}.pkl')
+    train_and_store.train_and_store(net=net, 
+            epoch_map=epoch_map,
+            number_of_epochs=number_of_epochs,
+            save_file=save_file)
+    
+
+if __name__ == '__main__':
+    for seed in seed_list:
+        # Tensorboard monitoring
+        writer = SummaryWriter(log_dir=f'/home/martin09/tmp/tensorboard/'\
+                f'run_eiv_yacht_lr_{lr:.4f}_seed'\
+                f'_{seed}_uregu_{unscaled_reg:.1f}_p_{p:.2f}'\
+                f'_fixed_std_x_{fixed_std_x:.3f}')
+        print(f'>>>>SEED: {seed}')
+        for init_std_y in init_std_y_list:
+            print(f'Using init_std_y={init_std_y:.3f}')
+            train_on_data(init_std_y, seed)
diff --git a/Experiments/train_noneiv_energy.py b/Experiments/train_noneiv_energy.py
index 04f299af8f0c92aeffd393f03379319844c248c8..106bc14e5802054f0bcfdbfdc3d7f4ebbd13d3ab 100644
--- a/Experiments/train_noneiv_energy.py
+++ b/Experiments/train_noneiv_energy.py
@@ -18,13 +18,13 @@ from EIVTrainingRoutines import train_and_store, loss_functions
 lr = 1e-3
 batch_size = 32
 test_batch_size = 600
-number_of_epochs = 300
+number_of_epochs = 600
 unscaled_reg = 10
 report_point = 5
 p = 0.2
-lr_update = 50
+lr_update = 100
 # pretraining = 300
-epoch_offset = 50
+epoch_offset = 100
 init_std_y_list = [0.5]
 gamma = 0.5
 hidden_layers = [1024, 1024, 1024, 1024]