diff --git a/examples/analytical-function/example_analytical_function.py b/examples/analytical-function/example_analytical_function.py
index 70124f9cf8b4681c02ec99141d0d5c2fc78df11d..30ff3fac4dfd3390135aed838295742d7fcef932 100644
--- a/examples/analytical-function/example_analytical_function.py
+++ b/examples/analytical-function/example_analytical_function.py
@@ -237,8 +237,9 @@ if __name__ == "__main__":
     PostPCE.check_accuracy(n_samples=300)
 
     # Compute the moments and compare with the Monte-Carlo reference
-    if MetaModelOpts.meta_model_type != 'GPE':
-        PostPCE.plot_moments()
+    # TODO: generalize the moment calculation
+    #if MetaModelOpts.meta_model_type != 'GPE':
+    #    PostPCE.plot_moments()
 
     # Plot the evolution of the KLD,BME, and Modified LOOCV error
     if MetaModelOpts.ExpDesign.method == 'sequential':
@@ -246,6 +247,7 @@ if __name__ == "__main__":
         PostPCE.plot_seq_design_diagnostics(refBME_KLD)
 
     # Plot the sobol indices
+    # TODO: generalize this by using moment calculation(?)
     if MetaModelOpts.meta_model_type != 'GPE':
         total_sobol = PostPCE.sobol_indices()
 
diff --git a/examples/user_guide/example_user_guide.py b/examples/user_guide/example_user_guide.py
index 4c191b7de270ef6a817bbf2fd3cc24b9d8d3f3a3..5a32c15d6d49493de221e02557e950036918b034 100644
--- a/examples/user_guide/example_user_guide.py
+++ b/examples/user_guide/example_user_guide.py
@@ -66,7 +66,7 @@ if __name__ == '__main__':
     MetaMod.pce_q_norm = 1
     
     # TODO: add check in Metamod/ExpDesign to ensure that n_init_samples matches the length of the given .X!
-    ExpDesign.n_init_samples = 30
+    ExpDesign.n_init_samples = 40
     ExpDesign.sampling_method = 'random'#'user'
     #ExpDesign.X = samples
     
diff --git a/src/bayesvalidrox/surrogate_models/engine.py b/src/bayesvalidrox/surrogate_models/engine.py
index 8d27310b70737a05f207314b30af0e29be280001..92e5dbf663973ba74cfbbea7bb7411f8c2605eb2 100644
--- a/src/bayesvalidrox/surrogate_models/engine.py
+++ b/src/bayesvalidrox/surrogate_models/engine.py
@@ -263,9 +263,13 @@ class Engine:
             pce = False
         #mc_ref = True if bool(self.Model.mc_reference) else False
         mc_ref = False
-        if self.Model.mc_reference != {}:
+        # TODO: somehow get issues from these comparisons, thus doubled it here
+        if self.Model.mc_reference is not {}:
             mc_ref = True
             self.Model.read_observation('mc_ref')
+        #elif self.Model.mc_reference != {}:
+        #    mc_ref = True
+        #    self.Model.read_observation('mc_ref')
 
         # Get the parameters
         max_n_samples = self.ExpDesign.n_max_samples
@@ -348,7 +352,7 @@ class Engine:
             for out_name in output_name:
                 y = self.ExpDesign.Y[out_name]
                 Scores_all.append(list(
-                    self.MetaModel.score_dict['b_1'][out_name].values()))
+                    self.MetaModel.LOOCV_score_dict['b_1'][out_name].values()))
                 if self.MetaModel.dim_red_method.lower() == 'pca':
                     pca = self.MetaModel.pca['b_1'][out_name]
                     components = pca.transform(y)
@@ -494,7 +498,7 @@ class Engine:
                         for out_name in output_name:
                             y = self.ExpDesign.Y[out_name]
                             Scores_all.append(list(
-                                self.MetaModel.score_dict['b_1'][out_name].values()))
+                                self.MetaModel.LOOCV_score_dict['b_1'][out_name].values()))
                             if self.MetaModel.dim_red_method.lower() == 'pca':
                                 pca = self.MetaModel.pca['b_1'][out_name]
                                 components = pca.transform(y)
@@ -1013,6 +1017,8 @@ class Engine:
         """
         if self.Model.mc_reference is {}:
             raise AttributeError('Model.mc_reference needs to be given to calculate the surrogate error!')
+        #elif self.Model.mc_reference == {}:
+        #    raise AttributeError('Model.mc_reference needs to be given to calculate the surrogate error!')
         
         # Compute the mean and std based on the MetaModel
         means, stds = self.MetaModel.compute_moments()
diff --git a/src/bayesvalidrox/surrogate_models/polynomial_chaos.py b/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
index fc483d57d4e324536ac9d69f29cfc4a9533bb166..103edfa8dae25b287bbd274e9a81be7c0aff2590 100644
--- a/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
+++ b/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
@@ -100,10 +100,10 @@ class PCE(MetaModel):
                  pce_reg_method='OLS', bootstrap_method='fast',
                  n_bootstrap_itrs=1, pce_deg=1, pce_q_norm=1.0,
                  dim_red_method='no', verbose=False):
-        # Check if the result would be gaussian
-        
+        # Check if the surrogate outputs gaussian results     
         is_gaussian = self.check_is_gaussian(pce_reg_method, n_bootstrap_itrs)
         
+        # Use parent init
         super().__init__(input_obj, bootstrap_method,
                      n_bootstrap_itrs, dim_red_method, is_gaussian, 
                      verbose)
@@ -121,22 +121,14 @@ class PCE(MetaModel):
         self._error_clf_poly = None
         self._errorRegMethod = None
         self._univ_p_val = None
-        #self.out_names = None
         self._allBasisIndices = None
         self._deg_array = None
-        #self.n_samples = None
-        #self.CollocationPoints = None # TODO: rename this!
-        #self.pca = None
         self._LCerror = None
         self._clf_poly = None
-        self._score_dict = None
         self._basis_dict = None
         self._coeffs_dict = None
         self._q_norm_dict = None
         self._deg_dict = None
-        #self.n_params = None
-        #self.ndim = None
-        #self.rmse = None
         
     def check_is_gaussian(self, _pce_reg_method, n_bootstrap_itrs) -> bool:
         """
@@ -202,7 +194,7 @@ class PCE(MetaModel):
         self._q_norm_dict = self.auto_vivification()
         self._coeffs_dict = self.auto_vivification()
         self._basis_dict = self.auto_vivification()
-        self._score_dict = self.auto_vivification()
+        self.LOOCV_score_dict = self.auto_vivification()
         self._clf_poly = self.auto_vivification()
         self._LCerror = self.auto_vivification()
         if self.dim_red_method.lower() == 'pca':
@@ -376,7 +368,7 @@ class PCE(MetaModel):
                     self._q_norm_dict[f'b_{b_i + 1}'][key][f"y_{i + 1}"] = out[i]['qnorm']
                     self._coeffs_dict[f'b_{b_i + 1}'][key][f"y_{i + 1}"] = out[i]['coeffs']
                     self._basis_dict[f'b_{b_i + 1}'][key][f"y_{i + 1}"] = out[i]['multi_indices']
-                    self._score_dict[f'b_{b_i + 1}'][key][f"y_{i + 1}"] = out[i]['LOOCVScore']
+                    self.LOOCV_score_dict[f'b_{b_i + 1}'][key][f"y_{i + 1}"] = out[i]['LOOCVScore']
                     self._clf_poly[f'b_{b_i + 1}'][key][f"y_{i + 1}"] = out[i]['_clf_poly']
                     # self._LCerror[f'b_{b_i+1}'][key][f"y_{i+1}"] = out[i]['_LCerror']
                     
@@ -567,8 +559,8 @@ class PCE(MetaModel):
         elif reg_method.lower() == 'ebl':
             _clf_poly = EBLinearRegression(optimizer='em')
 
-        print(X.shape)
-        print(y.shape)
+        #print(X.shape)
+        #print(y.shape)
         _clf_poly.fit(X, y)
 
         # Select the nonzero entries of coefficients
@@ -576,6 +568,7 @@ class PCE(MetaModel):
             nnz_idx = np.nonzero(_clf_poly.coef_)[0]
         else:
             nnz_idx = np.arange(_clf_poly.coef_.shape[0])
+        #print(nnz_idx)
 
         # This is for the case where all outputs are zero, thereby
         # all coefficients are zero
@@ -660,6 +653,8 @@ class PCE(MetaModel):
 
                 # Assemble the Psi matrix
                 Psi = create_psi(BasisIndices, self._univ_p_val)
+                #print(Psi.shape)
+                
 
                 # Calulate the cofficients of the metamodel
                 outs = self.regression(Psi, ED_Y, BasisIndices)
@@ -817,14 +812,12 @@ class PCE(MetaModel):
             method='user'
         )
         # Compute univariate bases for the given samples
-        _univ_p_val = None
         _univ_p_val = self.univ_basis_vals(
             samples,
             n_max=np.max(self._pce_deg)
         )
+        #print(_univ_p_val)
 
-        mean_pred = None
-        std_pred = None
         mean_pred_b = {}
         std_pred_b = {}
         b_i = 0
@@ -847,6 +840,7 @@ class PCE(MetaModel):
                     # Assemble Psi matrix
                     basis = self._basis_dict[f'b_{b_i + 1}'][output][in_key]
                     psi = create_psi(basis, _univ_p_val)
+                    #print(basis)
                     
                     #print(psi.shape)
 
@@ -854,6 +848,10 @@ class PCE(MetaModel):
                     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]
+                        #import inspect
+                        #args = inspect.getfullargspec(_clf_poly.predict)
+                        # TODO: check if 'return_std' in this 
+                        #if args.varargs is not None and 'return_std' in args.varargs:
                         try:
                             y_mean, y_std = _clf_poly.predict(
                                 psi, return_std=True
@@ -861,6 +859,9 @@ class PCE(MetaModel):
                         except TypeError:
                             y_mean = _clf_poly.predict(psi)
                             y_std = np.zeros_like(y_mean)
+                        #else:
+                        #    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]
diff --git a/src/bayesvalidrox/surrogate_models/sequential_design.py b/src/bayesvalidrox/surrogate_models/sequential_design.py
index a373f8dd16a716a670478ad2b1bf044619327290..38c5c3ac9df71c573bac6114aaa61e169e8718dc 100644
--- a/src/bayesvalidrox/surrogate_models/sequential_design.py
+++ b/src/bayesvalidrox/surrogate_models/sequential_design.py
@@ -1125,7 +1125,7 @@ class SequentialDesign:
         # TODO: Only one psi can be selected.
         # Suggestion: Go for the one with the highest LOO error
         # TODO: this is just a patch, need to look at again!
-        Scores = list(self.MetaModel.score_dict['b_1'][OutputName].values())
+        Scores = list(self.MetaModel.LOOCV_score_dict['b_1'][OutputName].values())
         ModifiedLOO = [1 - score for score in Scores]
         outIdx = np.argmax(ModifiedLOO)
 
@@ -1535,7 +1535,7 @@ class SequentialDesign:
 
         """
         # Compute the mean and std based on the MetaModel
-        pce_means, pce_stds = self.MetaModel._compute_pce_moments()
+        pce_means, pce_stds = self.MetaModel.compute_moments()
 
         # Compute the root mean squared error
         for output in self.out_names:
diff --git a/src/bayesvalidrox/surrogate_models/surrogate_models.py b/src/bayesvalidrox/surrogate_models/surrogate_models.py
index a8537f66692a80cbc26915839be271663e28468f..49733129f35c6c8d601b3ec00820094f1fbc3ac4 100644
--- a/src/bayesvalidrox/surrogate_models/surrogate_models.py
+++ b/src/bayesvalidrox/surrogate_models/surrogate_models.py
@@ -94,12 +94,11 @@ class MetaModel:
         self.CollocationPoints = None # TODO: rename this
         self.pca = None
         self.LCerror = None
-        #self.score_dict = None
-        #self.basis_dict = None
+        self.LOOCV_score_dict = None # TODO: generic function that calculates this score!?
         #self.x_scaler = None
         #self.gp_poly = None
         self.n_params = None
-        self.ndim = None
+        self.ndim = None # TODO: this should be the same as n_params!?
         #self.init_type = None # TODO: this parameter seems to not be used
         #self.rmse = None # TODO: check if this is needed