diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index 1748f30907b9795ef4bbc586838d6d370b5a0e0b..636688190f828d3cb327f5b19571d8c1456b6bae 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -866,10 +866,10 @@ class Metamodel:
     def PCATransformation(self, Output):
         
         # Transform via Principal Component Analysis
-        from sklearn.decomposition import PCA as sklearnPCA
+        from sklearn.decomposition import IncrementalPCA as sklearnPCA
         # from sklearn.decomposition import TruncatedSVD as TruncatedSVD
         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:
@@ -879,7 +879,7 @@ class Metamodel:
         
         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)
 
         return pca, OutputMatrix
@@ -893,10 +893,10 @@ class Metamodel:
                                               ConstantKernel)
 
 
-        kernels = [1.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-5, 1e5)),
-                   1.0 * RationalQuadratic(length_scale=0.2, alpha=1.0),
+        kernels = [np.var(y) * RBF(length_scale=1.0, length_scale_bounds=(1e-5, 1e5)),
+                   np.var(y) * RationalQuadratic(length_scale=0.2, alpha=1.0),
                    # 1.0 * ExpSineSquared(length_scale=1.0, length_scale_bounds=(1e-05, 1e05)),
-                   1.0 * Matern(length_scale=1.0, length_scale_bounds=(1e-5, 1e5),
+                   np.var(y) * Matern(length_scale=1.0, length_scale_bounds=(1e-5, 1e5),
                                 nu=1.5)]
         
         if autoSelect:# Automatic selection of the kernel
@@ -914,7 +914,7 @@ class Metamodel:
             gp = gp[np.argmax(BME)]    
         
         else:
-            gp = GaussianProcessRegressor(kernel=kernels[0], n_restarts_optimizer=3,
+            gp = GaussianProcessRegressor(kernel=kernels[-1], n_restarts_optimizer=3,
                                           normalize_y=True)
             gp.fit(X, y)
             
@@ -2186,7 +2186,7 @@ class Metamodel:
             
             # TODO: Plot PCE vs Components
             # PCA = self.pca[Outkey]
-            # self.plot_components(PCA, Y, PCEOutputs_mean, PCEOutputs_std, plotname='PCA_PCE_'+Outkey)
+            # self.plot_components(PCA, Y[Outkey], PCEOutputs_mean, PCEOutputs_std, plotname='PCA_PCE_'+Outkey)
                 
             if self.index is None:
                 if self.DimRedMethod.lower() == 'pca':