diff --git a/src/bayesvalidrox/bayes_inference/bayes_inference.py b/src/bayesvalidrox/bayes_inference/bayes_inference.py
index 5bc144960907c7dbdd8fffa2587453d9c8f246f3..7284348266fa88dd93a7380633437cc713e19291 100644
--- a/src/bayesvalidrox/bayes_inference/bayes_inference.py
+++ b/src/bayesvalidrox/bayes_inference/bayes_inference.py
@@ -302,7 +302,7 @@ class BayesInference:
                 samples = self.prior_samples
 
             # Take care of an additional Sigma2s
-            self.prior_samples = samples[:, :self.engine.MetaModel.n_params]
+            self.prior_samples = samples[:, :self.engine.MetaModel.ndim]
 
             # Update number of samples
             self.n_prior_samples = self.prior_samples.shape[0]
@@ -585,7 +585,7 @@ class BayesInference:
         # Evaluate the MetaModel
         if self.emulator:
             print('Evaluating the metamodel on the prior samples.')
-            y_hat, y_std = MetaModel.eval_metamodel(samples=self.prior_samples)
+            y_hat, y_std = MetaModel.eval_metamodel(self.prior_samples)
             self.__mean_pce_prior_pred = y_hat
             self._std_pce_prior_pred = y_std
 
diff --git a/src/bayesvalidrox/post_processing/post_processing.py b/src/bayesvalidrox/post_processing/post_processing.py
index 361a2c5f666d043f50e54b94ddd4441763d08176..0d7fbad16b8d25482f1d350466cb8e713373d157 100644
--- a/src/bayesvalidrox/post_processing/post_processing.py
+++ b/src/bayesvalidrox/post_processing/post_processing.py
@@ -252,7 +252,7 @@ class PostProcessing:
             outputs = self._eval_model(samples, key_str='validSet')
 
         # Run the PCE model with the generated samples
-        pce_outputs, _ = MetaModel.eval_metamodel(samples=samples)
+        pce_outputs, _ = MetaModel.eval_metamodel(samples)
 
         self.rmse = {}
         self.valid_error = {}
diff --git a/src/bayesvalidrox/surrogate_models/polynomial_chaos.py b/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
index 860eb80f5cc6ea1c2e3670df07c222f819b0b464..92aa33cc5d606d4a0a35ce9f978fa8487a2ddeb9 100644
--- a/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
+++ b/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
@@ -31,7 +31,7 @@ from .orthogonal_matching_pursuit import OrthogonalMatchingPursuit
 from .reg_fast_ard import RegressionFastARD
 from .reg_fast_laplace import RegressionFastLaplace
 
-from .surrogate_models import MetaModel, _preprocessing_fit, _bootstrap_fit
+from .surrogate_models import MetaModel, _preprocessing_fit, _bootstrap_fit, _preprocessing_eval, _bootstrap_eval
 
 from .supplementary import corr_loocv_error, create_psi
 
@@ -699,8 +699,12 @@ class PCE(MetaModel):
 
         return returnVars
 
+
+
     # -------------------------------------------------------------------------
-    def eval_metamodel(self, samples):
+    @_preprocessing_eval
+    @_bootstrap_eval
+    def eval_metamodel(self, samples, b_i = 0):
         """
         Evaluates metamodel at the requested samples. One can also generate
         nsamples.
@@ -717,96 +721,54 @@ class PCE(MetaModel):
         std_pred : dict
             Standard deviatioon of the predictions.
         """
-        # Transform into np array - can also be given as list
-        samples = np.array(samples)
-
-        # Transform samples to the independent space
-        samples = self.InputSpace.transform(samples, method='user')
-
         # Compute univariate bases for the given samples
         _univ_p_val = self.univ_basis_vals(samples,n_max=np.max(self._pce_deg))
 
-        mean_pred_b = {}
-        std_pred_b = {}
-        b_i = 0
-        # Loop over bootstrap iterations
-        for b_i in range(self.n_bootstrap_itrs):
-
-            # Extract model dictionary
-            model_dict = self._coeffs_dict[f'b_{b_i + 1}']
-
-            # Loop over outputs
-            mean_pred = {}
-            std_pred = {}
-            for output, values in model_dict.items():
-
-                mean = np.empty((len(samples), len(values)))
-                std = np.empty((len(samples), len(values)))
-                idx = 0
-                for in_key, _ in values.items():
-
-                    # Assemble Psi matrix
-                    basis = self._basis_dict[f'b_{b_i + 1}'][output][in_key]
-                    psi = create_psi(basis, _univ_p_val)
-
-                    # Prediction
-                    if self.bootstrap_method != 'fast' or b_i == 0:
-                        # with error bar, i.e. use _clf_poly
-                        _clf_poly = self._clf_poly[f'b_{b_i + 1}'][output][in_key]
-                        try:
-                            y_mean, y_std = _clf_poly.predict(
-                                psi, return_std=True
-                            )
-                        except TypeError:
-                            y_mean = _clf_poly.predict(psi)
-                            y_std = np.zeros_like(y_mean)
-                    else:
-                        # without error bar
-                        coeffs = self._coeffs_dict[f'b_{b_i + 1}'][output][in_key]
-                        y_mean = np.dot(psi, coeffs)
+        # Extract model dictionary
+        model_dict = self._coeffs_dict[f'b_{b_i + 1}']
+
+        # Loop over outputs
+        mean_pred = {}
+        std_pred = {}
+        for output, values in model_dict.items():
+            print(output, values)
+
+            mean = np.empty((len(samples), len(values)))
+            std = np.empty((len(samples), len(values)))
+            idx = 0
+            for in_key, _ in values.items():
+
+                # Assemble Psi matrix
+                basis = self._basis_dict[f'b_{b_i + 1}'][output][in_key]
+                psi = create_psi(basis, _univ_p_val)
+
+                # Prediction
+                if self.bootstrap_method != 'fast' or b_i == 0:
+                    # with error bar, i.e. use _clf_poly
+                    _clf_poly = self._clf_poly[f'b_{b_i + 1}'][output][in_key]
+                    try:
+                        y_mean, y_std = _clf_poly.predict(
+                            psi, return_std=True
+                        )
+                    except TypeError:
+                        y_mean = _clf_poly.predict(psi)
                         y_std = np.zeros_like(y_mean)
-
-                    mean[:, idx] = y_mean
-                    std[:, idx] = y_std
-                    idx += 1
-
-                # Save predictions for each output
-                if self.dim_red_method.lower() == 'pca':
-                    pca = self.pca[f'b_{b_i + 1}'][output]
-                    mean_pred[output] = pca.inverse_transform(mean)
-                    std_pred[output] = np.zeros(mean.shape)
                 else:
-                    mean_pred[output] = mean
-                    std_pred[output] = std
+                    # without error bar
+                    coeffs = self._coeffs_dict[f'b_{b_i + 1}'][output][in_key]
+                    y_mean = np.dot(psi, coeffs)
+                    y_std = np.zeros_like(y_mean)
 
-            # Save predictions for each bootstrap iteration
-            mean_pred_b[b_i] = mean_pred
-            std_pred_b[b_i] = std_pred
+                mean[:, idx] = y_mean
+                std[:, idx] = y_std
+                idx += 1
 
-        # Change the order of nesting
-        mean_pred_all = {}
-        for i in sorted(mean_pred_b):
-            for k, v in mean_pred_b[i].items():
-                if k not in mean_pred_all:
-                    mean_pred_all[k] = [None] * len(mean_pred_b)
-                mean_pred_all[k][i] = v
-
-        # Compute the moments of predictions over the predictions
-        for output in self.out_names:
-            # Only use bootstraps with finite values
-            finite_rows = np.isfinite(
-                mean_pred_all[output]).all(axis=2).all(axis=1)
-            outs = np.asarray(mean_pred_all[output])[finite_rows]
-
-            # Compute mean and stdev
-            mean_pred[output] = np.mean(outs, axis=0)
-            if self.n_bootstrap_itrs > 1:
-                std_pred[output] = np.std(outs, axis=0)
-            else:
-                std_pred[output] = std_pred_b[self.n_bootstrap_itrs-1][output]
+            mean_pred[output] = mean
+            std_pred[output] = std
 
         return mean_pred, std_pred
 
+
     # -------------------------------------------------------------------------
     def __select_degree(self):#, ndim, n_samples):
         """
diff --git a/src/bayesvalidrox/surrogate_models/surrogate_models.py b/src/bayesvalidrox/surrogate_models/surrogate_models.py
index 133c420ad82787e5f424f0d613eaef620e63dcbb..ddd57f89fe7bd26f28a7ac2eebcca9c33a5f233a 100644
--- a/src/bayesvalidrox/surrogate_models/surrogate_models.py
+++ b/src/bayesvalidrox/surrogate_models/surrogate_models.py
@@ -240,6 +240,106 @@ def transform_y(self, y, b_i=0, trafo_type = ''):
 
     return y_transform
 
+def _preprocessing_eval(eval_function):
+    """
+    This decorator performs the pre- and postprocessing for evaluating the
+    metamodel.
+
+    Parameters
+    ----------
+    eval_function : function
+        Function with the same signature as self.eval_metamodel().
+
+    Returns
+    -------
+    decorator : function
+        The decorated function.
+
+    """
+    @wraps(eval_function)
+    def decorator(self , *args, **kwargs) :
+        # Transform into np array - can also be given as list
+        samples = args[0]
+        samples = np.array(samples)
+
+        # Transform samples to the independent space
+        samples = self.InputSpace.transform(samples, method='user')
+
+        
+        out = eval_function(self, samples, **kwargs)
+        
+        return out
+    
+    return decorator
+
+
+def _bootstrap_eval(eval_function):
+    """
+    Decorator that applies bootstrap iterations around the evaluation of a 
+    MetaModel.
+    Also performs additional transformations to the training outputs, e.g. PCA.
+
+    Parameters
+    ----------
+    eval_function : function
+        Function with the same signature as self.eval_metamodel().
+
+    Returns
+    -------
+    decorator : function
+        The decorated function.
+
+    """
+    @wraps(eval_function)
+    def decorator(self , *args, **kwargs) :
+        mean_pred_b = {}
+        std_pred_b = {}
+        # Loop over bootstrap iterations
+        for b_i in range(self.n_bootstrap_itrs):
+            kwargs['b_i'] = b_i
+            mean_pred, std_pred = eval_function(self, *args, **kwargs)
+
+            # Apply inverse transformations
+            for i in range(mean_pred[list(mean_pred.keys())[0]].shape[1]):
+                # Save predictions for each output
+                if self.dim_red_method.lower() == 'pca':
+                    pca = self.pca[f'b_{b_i + 1}'][f"y_{i + 1}"]
+                    mean_pred[f"y_{i + 1}"] = pca.inverse_transform(mean_pred)
+                    std_pred[f"y_{i + 1}"] = np.zeros(mean_pred.shape)
+                else:
+                    mean_pred[f"y_{i + 1}"] = mean_pred
+                    std_pred[f"y_{i + 1}"] = std_pred
+
+            # Save predictions for each bootstrap iteration
+            mean_pred_b[b_i] = mean_pred
+            std_pred_b[b_i] = std_pred
+
+        # Change the order of nesting
+        mean_pred_all = {}
+        for i in sorted(mean_pred_b):
+            for k, v in mean_pred_b[i].items():
+                if k not in mean_pred_all:
+                    mean_pred_all[k] = [None] * len(mean_pred_b)
+                mean_pred_all[k][i] = v
+
+        # Compute the moments of predictions over the predictions
+        for output in self.out_names:
+            # Only use bootstraps with finite values
+            finite_rows = np.isfinite(
+                mean_pred_all[output]).all(axis=2).all(axis=1)
+            outs = np.asarray(mean_pred_all[output])[finite_rows]
+
+            # Compute mean and stdev
+            mean_pred[output] = np.mean(outs, axis=0)
+            if self.n_bootstrap_itrs > 1:
+                std_pred[output] = np.std(outs, axis=0)
+            else:
+                std_pred[output] = std_pred_b[self.n_bootstrap_itrs-1][output]
+
+        return mean_pred, std_pred
+    
+    return decorator
+
 
 class MetaModel:
     """