diff --git a/src/bayesvalidrox/bayes_inference/bayes_model_comparison.py b/src/bayesvalidrox/bayes_inference/bayes_model_comparison.py
index a10e61f555dc4be4433224ed802b9ed4206cd521..21a0fe30299ed92fb305418c50095ce0336833e5 100644
--- a/src/bayesvalidrox/bayes_inference/bayes_model_comparison.py
+++ b/src/bayesvalidrox/bayes_inference/bayes_model_comparison.py
@@ -20,10 +20,11 @@ plt.style.use(os.path.join(os.path.split(__file__)[0],
 
 
 class BayesModelComparison:
-    def __init__(self, justifiability=True, n_bootstarp=1000,
-                 data_noise_level=0.01, just_n_meas=2):
+    def __init__(self, justifiability=True, perturbed_data=None,
+                 n_bootstarp=1000, data_noise_level=0.01, just_n_meas=2):
 
         self.justifiability = justifiability
+        self.perturbed_data = perturbed_data
         self.n_bootstarp = n_bootstarp
         self.data_noise_level = data_noise_level
         self.just_n_meas = just_n_meas
@@ -143,17 +144,18 @@ class BayesModelComparison:
         out_names = metaModel.ModelObj.Output.names
 
         # Perturb observations for Bayes Factor
-        perturbed_data = self.__perturb_data(
-            metaModel.ModelObj.observations, out_names, n_bootstarp,
-            noise_level=self.data_noise_level)
+        if self.perturbed_data is None:
+            self.perturbed_data = self.__perturb_data(
+                    metaModel.ModelObj.observations, out_names, n_bootstarp,
+                    noise_level=self.data_noise_level)
 
         # Only for Bayes Factor
         if not justifiability:
-            return perturbed_data
+            return self.perturbed_data
 
         # Generate data
         for i in range(n_bootstarp):
-            y_data = perturbed_data[i].reshape(1, -1)
+            y_data = self.perturbed_data[i].reshape(1, -1)
             justData = np.tril(np.repeat(y_data, y_data.shape[1], axis=0))
             # Use surrogate runs for data-generating process
             for key, metaModel in modelDict.items():