diff --git a/.gitignore b/.gitignore
index acb906e2023293020272c9329370d5d02129c7ce..7aa5fc0a2f338b13b8bd7f9765e9a051d98a96a0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -65,3 +65,54 @@ BayesValidRox/tests/AnalyticalFunction/Outputs_sequential_BayesOptDesign/
 BayesValidRox/tests/AnalyticalFunction/Outputs_sequential_alphabetic/
 BayesValidRox/tests/AnalyticalFunction/Polynomials/
 BayesValidRox/tests/AnalyticalFunction/__pycache__/
+BayesValidRox/BayesInference/__pycache__/BayesInference.cpython-37.pyc
+BayesValidRox/BayesInference/__pycache__/Discrepancy.cpython-37.pyc
+BayesValidRox/BayesInference/__pycache__/__init__.cpython-37.pyc
+BayesValidRox/PostProcessing/__pycache__/PostProcessing.cpython-37.pyc
+BayesValidRox/PostProcessing/__pycache__/__init__.cpython-37.pyc
+BayesValidRox/PyLink/__pycache__/FuncForwardModel.cpython-37.pyc
+BayesValidRox/PyLink/__pycache__/PyLinkForwardModel.cpython-37.pyc
+BayesValidRox/PyLink/__pycache__/__init__.cpython-37.pyc
+BayesValidRox/surrogate_models/__pycache__/ExpDesigns.cpython-37.pyc
+BayesValidRox/surrogate_models/__pycache__/Exploration.cpython-37.pyc
+BayesValidRox/surrogate_models/__pycache__/Input.cpython-37.pyc
+BayesValidRox/surrogate_models/__pycache__/Voronoi.cpython-37.pyc
+BayesValidRox/surrogate_models/__pycache__/__init__.cpython-37.pyc
+BayesValidRox/surrogate_models/__pycache__/aPoly_Construction.cpython-37.pyc
+BayesValidRox/surrogate_models/__pycache__/surrogate_models.cpython-37.pyc
+BayesValidRox/tests/AnalyticalFunction/Outputs_demo/
+BayesValidRox/tests/BeamTest/Beam_Results.pkl
+BayesValidRox/tests/BeamTest/Polynomials/Roots_degree7_Pa1.npy
+BayesValidRox/tests/BeamTest/Polynomials/Roots_degree7_Pa2.npy
+BayesValidRox/tests/BeamTest/Polynomials/Roots_degree7_Pa3.npy
+BayesValidRox/tests/BeamTest/Polynomials/Roots_degree7_Pa4.npy
+BayesValidRox/tests/BeamTest/__pycache__/read_Beam_Deflection.cpython-37.pyc
+BayesValidRox/tests/BeamTest/surrogateBeam/Outputs_PostPorecessing/
+BayesValidRox/tests/BeamTest/surrogateBeam/Polynomials/
+BayesValidRox/tests/BeamTest/surrogateBeam/__pycache__/read_Beam_Deflection.cpython-37.pyc
+BayesValidRox/tests/BraninTest/Branin_Results.pkl
+BayesValidRox/tests/BraninTest/Outputs_Bayes_Calib/
+BayesValidRox/tests/BraninTest/Outputs_PostPorecessing/
+BayesValidRox/tests/BraninTest/Outputs_SeqPosteriorComparison/
+BayesValidRox/tests/BraninTest/Outputs_normal_BayesOptDesign/
+BayesValidRox/tests/BraninTest/Outputs_sequential_BayesOptDesign/
+BayesValidRox/tests/BraninTest/Outputs_sequential_dual annealing/
+BayesValidRox/tests/BraninTest/__pycache__/Branin.cpython-37.pyc
+BayesValidRox/tests/Co2BenchmarkTest/CollocationPoints.npy
+BayesValidRox/tests/Co2BenchmarkTest/ExpDesign_Y.npy
+BayesValidRox/tests/Co2BenchmarkTest/Results/
+"BayesValidRox/tests/Co2BenchmarkTest/uqBenchmarkData/K\303\266ppel2019_Article_ComparisonOfData-drivenUncerta.pdf"
+BayesValidRox/tests/Co2BenchmarkTest/uqBenchmarkData/mc_10.zip
+BayesValidRox/tests/IshigamiTest/Outputs_PostPorecessing/
+BayesValidRox/tests/IshigamiTest/Polynomials/
+BayesValidRox/tests/MECBM/.model_exe_orig.sh.swp
+BayesValidRox/tests/MECBM/SurrogateModel/MECBM.zip
+BayesValidRox/tests/MECBM/SurrogateModel/MECBMValid.zip
+BayesValidRox/tests/MECBM/SurrogateModel/MECBM_SA_Results.pkl
+BayesValidRox/tests/MECBM/SurrogateModel/NormalDesign_100/
+BayesValidRox/tests/MECBM/SurrogateModel/SeqDesign_25_100/
+BayesValidRox/tests/MECBM/SurrogateModel/SeqDesign_25_100_APP/
+BayesValidRox/tests/MECBM/SurrogateModel/SeqDesign_25_100_DPP/
+BayesValidRox/tests/MECBM/SurrogateModel/SeqDesign_25_150_DKL/
+BayesValidRox/tests/MECBM/SurrogateModel/__pycache__/read_MECBM.cpython-37.pyc
+BayesValidRox/tests/MECBM/model_exe_orig.sh
diff --git a/BayesValidRox/PostProcessing/PostProcessing.py b/BayesValidRox/PostProcessing/PostProcessing.py
index 660f8e1a2cb332d43d51f9d71a0d265fb4f1952c..c6f9dde320f468e7f1c9533a2339091407e3f04e 100644
--- a/BayesValidRox/PostProcessing/PostProcessing.py
+++ b/BayesValidRox/PostProcessing/PostProcessing.py
@@ -54,7 +54,7 @@ class PostProcessing:
                 # Mean = c_0
                 PCEMean[idx] = coeffs[0]
                 # Std = sum(coeffs[1:]^2)
-                PCEStd[idx] = np.sum(np.square(coeffs[1:]))    
+                PCEStd[idx] = np.sqrt(np.sum(np.square(coeffs[1:])))
                 
             self.PCEMeans[Outkey] = PCEMean
             self.PCEStd[Outkey] = PCEStd
@@ -542,33 +542,33 @@ class PostProcessing:
         """
         PCEModel = self.PCEModel
         
-        plotList = ['Modified LOO error', 'KLD', 'BME']
-        seqList = [PCEModel.SeqModifiedLOO, PCEModel.SeqKLD, PCEModel.SeqBME]
-        
-        # Plot the evolution of the diagnostic criteria of the Sequential Experimental Design.
-        
-        markers = ['x', 'o', 'd']
-        colors = ['k', 'r', 'b']
+        plotList = ['Modified LOO error', 'KLD', 'BME', 'RMSEMean', 'RMSEStd']
+        seqList = [PCEModel.SeqModifiedLOO, PCEModel.SeqKLD, PCEModel.SeqBME, 
+                   PCEModel.seqRMSEMean, PCEModel.seqRMSEStd]
         
+        markers = ('x', 'o', 'd', '*', '+')
+        colors = ('k', 'darkgreen', 'b', 'navy', 'darkred')
         
+        # Plot the evolution of the diagnostic criteria of the Sequential Experimental Design.
         for plotidx, plot in enumerate(plotList): 
-            fig = plt.figure(figsize=(24, 16))
             
             SeqDict = seqList[plotidx]
             nameUtil = list(SeqDict.keys())
             
+            if len(nameUtil) == 0:
+                continue
+            
+            fig = plt.figure(figsize=(24, 16))
             
             for idx, name in enumerate(nameUtil):
                 SeqValues = SeqDict[name]
                 x_idx = np.arange(NrofInitSamples, SeqValues.shape[0]+NrofInitSamples)
-                marker = markers[idx]
-                color = colors[idx]
-                
-                if plot == 'Modified LOO error':
-                    plt.semilogy(x_idx, SeqValues, marker=marker, color=color, label=name)
-                else:
-                    plt.plot(x_idx, SeqValues, marker=marker, color=color, label=name)
 
+                if plot == 'KLD' or plot == 'BME':
+                    plt.plot(x_idx, SeqValues, marker=markers[idx], color=colors[idx], label=name)
+                else:
+                    plt.semilogy(x_idx, SeqValues, marker=markers[idx], color=colors[idx], label=name)
+                    
             plt.xlabel('Number of runs')
             plt.ylabel(plot)
             plt.legend()
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index ba5044ccf9d6a310e1ac8e99c9273a68cc0cea4c..1f74984535f141d4593dcf57e15bc3e04ab61f54 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -69,6 +69,8 @@ class aPCE:
         self.SeqModifiedLOO = []
         self.SeqBME = {}
         self.SeqKLD = {}
+        self.seqRMSEMean = {}
+        self.seqRMSEStd = {}
         self.Observations = []
         
         self.Likelihoods = None
@@ -1212,11 +1214,11 @@ class aPCE:
                 # Mean = c_0
                 PCEMean[idx] = coeffs[0]
                 # Std = sum(coeffs[1:]^2)
-                PCEStd[idx] = np.sum(np.square(coeffs[1:]))    
+                PCEStd[idx] = np.sqrt(np.sum(np.square(coeffs[1:])))
             
             # Compute the error between mean and std of PCEModel and OrigModel
-            MSE_Mean = mean_squared_error(df_MCReference['mean'], PCEMean)
-            MSE_std = mean_squared_error(df_MCReference['std'], PCEStd)
+            MSE_Mean = np.sqrt(mean_squared_error(df_MCReference['mean'], PCEMean))
+            MSE_std = np.sqrt(mean_squared_error(df_MCReference['std'], PCEStd))
             
             PCEMeans[Outkey] = PCEMean
             PCEStds[Outkey] = PCEStd
@@ -1283,7 +1285,10 @@ class aPCE:
             print('Posterior snapshot (initial) is being plotted...')
             self.posteriorPlot(initPosterior, MAP, parNames, 'SeqPosterior_init')
                    
-        
+        # Check the convergence of the Mean&Std
+        if self.ModelObj.MCReference:
+                initRMSEMean, initRMSEStd = self.error_Mean_Std()
+                print("Initial Mean and Std error: %s, %s"%(initRMSEMean, initRMSEStd))
         
         # Read the initial experimental design
         from copy import deepcopy
@@ -1300,6 +1305,8 @@ class aPCE:
         
         self.SeqModifiedLOO = {}
         self.SeqBME = {}
+        self.seqRMSEMean = {}
+        self.seqRMSEStd = {}
         
         for UtilityFunction in UtilityFunctions:
             
@@ -1320,6 +1327,9 @@ class aPCE:
             SeqModifiedLOO = np.array(ModifiedLOO)
             SeqBME = np.array([initBME])
             SeqKLD = np.array([initKLD])
+            if self.ModelObj.MCReference:
+                seqRMSEMean = np.array([initRMSEMean])
+                seqRMSEStd = np.array([initRMSEStd])
             
             postcnt = 1
             
@@ -1395,20 +1405,33 @@ class aPCE:
                     
                 postcnt += 1
                 
+                # Check the convergence of the Mean&Std
+                if self.ModelObj.MCReference:
+                    print('\n')
+                    RMSE_Mean, RMSE_std = self.error_Mean_Std()
+                    print("Updated Mean and Std error: %s, %s"%(RMSE_Mean, RMSE_std))
+                    print('\n')
+                    
                 print('\n')
                 print("Updated BME:", BME)
                 print("Updated KLD:", KLD)
                 print('\n')
                 
-                # Store the updated BME & KLD
+                # Store the updated BME & KLD & 
                 SeqBME  = np.vstack((SeqBME, BME))
                 SeqKLD  = np.vstack((SeqKLD, KLD))
+                if self.ModelObj.MCReference:
+                    seqRMSEMean = np.vstack((seqRMSEMean, RMSE_Mean))
+                    seqRMSEStd = np.vstack((seqRMSEStd, RMSE_std))
                 
             # Store updated ModifiedLOO and BME in dictonary
             self.SeqModifiedLOO[UtilityFunction] = SeqModifiedLOO
             self.SeqBME[UtilityFunction] = SeqBME
             self.SeqKLD[UtilityFunction] = SeqKLD
-            
+            if self.ModelObj.MCReference:
+                self.seqRMSEMean[UtilityFunction] = seqRMSEMean
+                self.seqRMSEStd[UtilityFunction] = seqRMSEStd
+                
         return PCEModel
     
     #--------------------------------------------------------------------------------------------------------
@@ -1432,7 +1455,6 @@ class aPCE:
         # TODO: Loop over outputs
         OutputName = Model.Output.Names
         
-        
         # Read the initial experimental design
         from copy import deepcopy
         Xinit = PCEModel.ExpDesign.X
@@ -1446,7 +1468,15 @@ class aPCE:
         Scores = [item for sublist in Scores_all for item in sublist]
         initModifiedLOO = [1-score for score in Scores]
         
+        # Check the initial EMSE error of the Mean&Std
+        if self.ModelObj.MCReference:
+                initRMSEMean, initRMSEStd = self.error_Mean_Std()
+                print("Initial Mean and Std error: %s, %s"%(initRMSEMean, initRMSEStd))
+                
         self.SeqModifiedLOO = {}
+        self.seqRMSEMean = {}
+        self.seqRMSEStd = {}
+        
         for UtilityFunction in UtilityFunctions:
             
             print ('>>>> UtilityFunction: %s<<<<'%UtilityFunction)
@@ -1463,11 +1493,10 @@ class aPCE:
             print("Initial ModifiedLOO:", initModifiedLOO)
             ModifiedLOO = initModifiedLOO
             SeqModifiedLOO = np.array(ModifiedLOO)
-            
-            # Check the convergence of the Mean&Std
             if self.ModelObj.MCReference:
-                MSE_Mean, MSE_std = self.error_Mean_Std()
-                print("Initial Mean and Std error: %s, %s"%(MSE_Mean, MSE_std))
+                seqRMSEMean = np.array([initRMSEMean])
+                seqRMSEStd = np.array([initRMSEStd])
+            
             
             # Start Sequential Experimental Design
             while (TotalNSamples < MaxNSamples) and any(LOO > ModifiedLOOThreshold for LOO in ModifiedLOO):
@@ -1530,13 +1559,18 @@ class aPCE:
                 # Check the convergence of the Mean&Std
                 if self.ModelObj.MCReference:
                     print('\n')
-                    MSE_Mean, MSE_std = self.error_Mean_Std()
-                    print("Updated Mean and Std error: %s, %s"%(MSE_Mean, MSE_std))
+                    RMSE_Mean, RMSE_std = self.error_Mean_Std()
+                    print("Updated Mean and Std error: %s, %s"%(RMSE_Mean, RMSE_std))
                     print('\n')
+                    seqRMSEMean = np.vstack((seqRMSEMean, RMSE_Mean))
+                    seqRMSEStd = np.vstack((seqRMSEStd, RMSE_std))
             
             # Store updated ModifiedLOO and BME in dictonary
             self.SeqModifiedLOO[UtilityFunction] = SeqModifiedLOO
-            
+            if self.ModelObj.MCReference:
+                self.seqRMSEMean[UtilityFunction] = seqRMSEMean
+                self.seqRMSEStd[UtilityFunction] = seqRMSEStd
+        
         return PCEModel
     
     #--------------------------------------------------------------------------------------------------------
diff --git a/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py b/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
index 400359bdbfe538d1a688df5a0a02dc94d3a0e183..03dbdb21312541f11123b9e977e1fac2f3535ce9 100644
--- a/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
+++ b/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
@@ -102,14 +102,14 @@ if __name__ == "__main__":
     # hypercube sampling of the input model or user-defined values of X and/or Y:
     MetaModelOpts.addExpDesign()
     
-    NrofInitSamples = 5
+    NrofInitSamples = 15
     MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
-    MetaModelOpts.ExpDesign.SamplingMethod = 'MC' # 1)MC 2)LHS 3)PCM 4)LSCM 5)user
+    MetaModelOpts.ExpDesign.SamplingMethod = 'LHS' # 1)MC 2)LHS 3)PCM 4)LSCM 5)user
     MetaModelOpts.ExpDesign.Method = 'sequential'  # 1) normal  2) sequential
     
     # Sequential experimental design (needed only for sequential ExpDesign)
     MetaModelOpts.ExpDesign.MaxNSamples = 50
-    MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-3
+    MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-5
     
     # Plot the posterior snapshots for SeqDesign
 #    MetaModelOpts.ExpDesign.PostSnapshot = True
diff --git a/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py b/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
index 49c525e0bdfc899dd05c3551c95d99d270946e6c..27f39d08efc0014b4dff5baa7dd1bccced6e08f1 100644
--- a/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
+++ b/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
@@ -203,13 +203,15 @@ if __name__ == "__main__":
     
     # ----- Define the discrepancy model -------
     obsData = pd.read_csv(Model.MeasurementFile)['CO2 saturation[-]']
+    error = 0.02*obsData.to_numpy().reshape(1,-1)[0]
+    error[0] = 1e-5
     # (Option A)
     # Do nothing
     
     # (Option B) # 0.2% of the data
     DiscrepancyOpts = Discrepancy('')
     DiscrepancyOpts.Type = 'Gaussian'
-    DiscrepancyOpts.Parameters = 0.02*obsData.to_numpy().reshape(1,-1)[0]     
+    DiscrepancyOpts.Parameters = error    
     BayesOpts.Discrepancy = DiscrepancyOpts
 
     # (Option C)