diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index 0d7fc81c20b6c43c7e098533ce1919fdd11e1c75..71827688944c4b1f4c4e63f47f8f6384e21caa32 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -158,7 +158,7 @@ class BayesInference:
         TotalOutputs = np.empty((SampleSize,0))
         for outputType in OutputType:
             TotalOutputs = np.hstack((TotalOutputs, Outputs[outputType]))
-        
+
         # Add the std of the PCE is chosen as emulator.
         if self.emulator:
             covMatrix_PCE = np.zeros((NofMeasurements, NofMeasurements), float)
@@ -178,7 +178,7 @@ class BayesInference:
             logL = multivariate_normal.logpdf(TotalOutputs, mean=Data, cov=covMatrix)
         except:
             logL = -np.inf
-
+            
         return logL
     
     #--------------------------------------------------------------------------------------------------------
diff --git a/BayesValidRox/PostProcessing/PostProcessing.py b/BayesValidRox/PostProcessing/PostProcessing.py
index 334c9ad7d375cfec98da22b278609ea0a8d51add..c50be9a86038dbeafb77f2687949a9926cfa84f6 100644
--- a/BayesValidRox/PostProcessing/PostProcessing.py
+++ b/BayesValidRox/PostProcessing/PostProcessing.py
@@ -67,7 +67,7 @@ class PostProcessing:
             
             if PCEModel.DimRedMethod.lower() == 'pca':
                 PCA = PCEModel.pca[Outkey]
-                self.PCEMeans[Outkey] = PCA.inverse_transform(PCEMean)
+                self.PCEMeans[Outkey] = PCA.mean_ + np.dot(PCEMean, PCA.components_)
                 self.PCEStd[Outkey] = np.sqrt(np.dot(PCEVar, PCA.components_**2))
             else:
                 self.PCEMeans[Outkey] = PCEMean
diff --git a/BayesValidRox/surrogate_models/ExpDesigns.py b/BayesValidRox/surrogate_models/ExpDesigns.py
index 1a1ed840f82bdefff4cab25ef26e3f263d03c2ef..ccb9d5d7ebd492e0f5c8142a2a3c162be9c3b609 100644
--- a/BayesValidRox/surrogate_models/ExpDesigns.py
+++ b/BayesValidRox/surrogate_models/ExpDesigns.py
@@ -31,7 +31,8 @@ class ExpDesigns:
         self.NCandidate = 1
         self.PostSnapshot = False
         self.stepSnapshot = 2
-        self.UtilityFunction = 'None'
+        self.ExploitMethod = 'Space-filling'
+        self.UtilityFunction = 'Space-filling'
         self.ExploreMethod = 'random'
         self.NrofCandGroups = 4
         self.nReprications = 1
diff --git a/BayesValidRox/surrogate_models/RegressionFastARD.py b/BayesValidRox/surrogate_models/RegressionFastARD.py
index 880baf488fce50a3c8209628f006df52c4d707b9..93819b75a89cc738ea4511550cc45662292ca909 100755
--- a/BayesValidRox/surrogate_models/RegressionFastARD.py
+++ b/BayesValidRox/surrogate_models/RegressionFastARD.py
@@ -236,7 +236,6 @@ class RegressionFastARD(LinearModel,RegressorMixin):
             # in case of almost perfect multicollinearity between some features
             # start from feature 0
             if np.sum( XXd - X_mean**2 < np.finfo(np.float32).eps ) > 0:#>= 0:
-                print("start from feature 0.")
                 A[0]       = np.finfo(np.float16).eps
                 active[0]  = True
             
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index 636688190f828d3cb327f5b19571d8c1456b6bae..791b6dff82798b13f9dcb093a200d60338d4b1b2 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -684,8 +684,9 @@ class Metamodel:
         LCerror = AllLCerror[str(bestIdx)]
         degree = DegreeArray[np.nanargmax(scores)]
         qnorm = float(qnorm[bestqIdx])
-
-        if verbose:
+        
+        # ------------------ Summary of results ------------------
+        if self.DisplayFlag:
             # Create PSI_Sparse by removing redundent terms
             nnz_idx = np.nonzero(coeffs)[0]
             BasisIndices_Sparse = PolynomialDegrees[nnz_idx]
@@ -697,8 +698,8 @@ class Metamodel:
                                                 round(len(BasisIndices_Sparse)/P, 3)))
             print('Final ModLOO error estimate: %.3e'%(1-max(scores)))
             print('\n'+'-'*50)
-        # ------------------ Summary of results ------------------
-        if self.DisplayFlag:
+        
+        if verbose:
             print('='*50)
             print(' '* 10 +' Summary of results ')
             print('='*50)
@@ -867,9 +868,10 @@ class Metamodel:
         
         # Transform via Principal Component Analysis
         from sklearn.decomposition import IncrementalPCA as sklearnPCA
-        # from sklearn.decomposition import TruncatedSVD as TruncatedSVD
+        # from sklearn.decomposition import TruncatedSVD as sklearnPCA
+
         n_samples, n_features = Output.shape
-        covar_matrix = sklearnPCA(n_components=None)#, svd_solver='full')
+        covar_matrix = sklearnPCA(n_components=None) #, svd_solver='full')
         covar_matrix.fit(Output)
         var = np.cumsum(np.round(covar_matrix.explained_variance_ratio_, decimals=5)*100)
         try:
@@ -878,14 +880,14 @@ class Metamodel:
             selected_n_components = min(n_samples, n_features)
         
         nComponents = min(n_samples, n_features, selected_n_components)
-
-        pca = sklearnPCA(n_components=nComponents)#, svd_solver='full')
+        pca = sklearnPCA(n_components=nComponents) #, svd_solver='full')
         OutputMatrix = pca.fit_transform(Output)
 
+        # pca.mean_ = np.zeros(Output.shape[1])
         return pca, OutputMatrix
 
     #--------------------------------------------------------------------------------------------------------
-    def GaussianProcessEmulator(self, X, y, autoSelect=True, varIdx=None):
+    def GaussianProcessEmulator(self, X, y, autoSelect=False, varIdx=None):
         
         from sklearn.gaussian_process import GaussianProcessRegressor
         from sklearn.gaussian_process.kernels import (RBF, Matern, RationalQuadratic,
@@ -914,7 +916,7 @@ class Metamodel:
             gp = gp[np.argmax(BME)]    
         
         else:
-            gp = GaussianProcessRegressor(kernel=kernels[-1], n_restarts_optimizer=3,
+            gp = GaussianProcessRegressor(kernel=kernels[1], n_restarts_optimizer=3,
                                           normalize_y=True)
             gp.fit(X, y)
             
@@ -1074,7 +1076,7 @@ class Metamodel:
         #                 prevBasisIndices = self.BasisDict[key]["y_"+str(idx)] if idx != 0 else None
                             
         #                 outDict = self.adaptiveSparseaPCE(CollocationPoints,OutputMatrix[:,idx],
-        #                                           idx,prevBasisIndices,verbose=True)
+        #                                           idx,prevBasisIndices)
                         
                         
                         
@@ -2183,7 +2185,7 @@ class Metamodel:
                 PCEOutputs_mean[:, idx] = y_mean
                 PCEOutputs_std[:, idx] = y_std
                 idx += 1
-            
+
             # TODO: Plot PCE vs Components
             # PCA = self.pca[Outkey]
             # self.plot_components(PCA, Y[Outkey], PCEOutputs_mean, PCEOutputs_std, plotname='PCA_PCE_'+Outkey)
@@ -2251,13 +2253,12 @@ class Metamodel:
         np.fill_diagonal(covMatrix_PCE, varPCE)
         
         # Aggregate the cov matrices
-        covMatrix = covMatrix + covMatrix_PCE
+        covMatrix += covMatrix_PCE
         
         # Compute likelihood
         self.Likelihoods = stats.multivariate_normal.pdf(TotalOutputs, 
                                                          mean=TotalData, 
-                                                         cov=covMatrix)#,
-                                                         # allow_singular=True)
+                                                         cov=covMatrix)
         
         return self.Likelihoods