diff --git a/BayesValidRox/PostProcessing/PostProcessing.py b/BayesValidRox/PostProcessing/PostProcessing.py
index 0f7993af0643a25e3a35b7eb273253a01c4678c2..8bd67e5a564556cc4bac07212c1fc05f06e4915f 100644
--- a/BayesValidRox/PostProcessing/PostProcessing.py
+++ b/BayesValidRox/PostProcessing/PostProcessing.py
@@ -11,6 +11,7 @@ import math
 import os
 from itertools import combinations, cycle
 import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
 plt.rcParams.update({'font.size': 24})
 plt.rc('figure', figsize = (24, 16))
 plt.rc('font', family='serif', serif='Arial')
@@ -19,8 +20,12 @@ plt.rc('text', usetex=True)
 plt.rc('xtick', labelsize=24)
 plt.rc('ytick', labelsize=24)
 plt.rc('axes', labelsize=24)
+plt.rc('axes', linewidth=2)
+plt.rc('axes', grid=True)
+plt.rc('grid', linestyle="-")
 plt.rc('savefig', dpi=1000)
 
+
 from sklearn.linear_model import LinearRegression
 from sklearn.metrics import mean_squared_error, r2_score
 
@@ -672,11 +677,12 @@ class PostProcessing:
                     if plot == 'KLD' or plot == 'BME':
                         # BME convergence if refBME is provided
                         if refBME_KLD is not None:
-                            if plot == 'BME': refValue, plotLabel = refBME_KLD[0], r'$BME-BME^{Ref.}$'
-                            if plot == 'KLD': refValue, plotLabel = refBME_KLD[1], r'$D_{KL}[p(\theta|y_*), p(\theta)] - D_{KL}^{Ref.}[p(\theta|y_*), p(\theta)]$'
+                            if plot == 'BME': refValue, plotLabel = refBME_KLD[0], r'$BME/BME^{Ref.}$'
+                            if plot == 'KLD': refValue, plotLabel = refBME_KLD[1], r'$D_{KL}[p(\theta|y_*),p(\theta)] / D_{KL}^{Ref.}[p(\theta|y_*), p(\theta)]$'
                             
                             # Difference between BME/KLD and the ref. values
-                            allErrors = allErrors - np.full((allErrors.shape), refValue)
+                            allErrors = np.divide(allErrors,np.full((allErrors.shape), refValue))
+                            
                             # Plot baseline for zero, i.e. no difference
                             plt.axhline(y=0.0, xmin=0, xmax=1, c='green', ls='--', lw=2)
                     
@@ -685,7 +691,12 @@ class PostProcessing:
                     draw_plot(allErrors[:,::step2], labels, edge_color, fill_color, idx)
                     
                  
-                plt.xticks(labels, labels)
+                # plt.xticks(labels, labels)
+                # Set the major and minor locators
+                ax.xaxis.set_major_locator(ticker.AutoLocator())
+                ax.xaxis.set_minor_locator(ticker.AutoMinorLocator())
+                ax.xaxis.grid(True, which='major', linestyle='-')
+                ax.xaxis.grid(True, which='minor', linestyle='--')
                 
                 from matplotlib.patches import Patch
                 legend_elements = [Patch(facecolor=fill_color[idx], edgecolor=edge_color[idx],
@@ -713,14 +724,14 @@ class PostProcessing:
                 
                 
             else:    
-                fig = plt.figure(figsize=(24, 16))
+                fig,ax = plt.subplots()
                 
                 for idx, name in enumerate(nameUtil):
                     SeqValues = SeqDict[name]
                     step = PCEModel.ExpDesign.NrofNewSample if PCEModel.ExpDesign.NrofNewSample!=1 else 1
                     x_idx = np.arange(NrofInitSamples, totalNSamples+1, step)
                     x_idx = np.hstack((x_idx, totalNSamples)) if totalNSamples not in x_idx else x_idx
-                    
+
                     if plot == 'KLD' or plot == 'BME':
                         # BME convergence if refBME is provided
                         if refBME_KLD is not None:
@@ -751,17 +762,19 @@ class PostProcessing:
                             plt.semilogy(x_idx, SeqValues[:,i], marker=markers[idx],
                                       ls='--', lw=2, color=colors[idx], alpha=0.15)
                         
-                        # plt.fill_between(x_idx, SeqValues.mean(axis=1)-2*np.std(SeqValues,axis=1), 
-                        #                  SeqValues.mean(axis=1)+2*np.std(SeqValues,axis=1), color=colors[idx],alpha=0.15)
-                
-                        
-                        plt.semilogy(x_idx, SeqValues.mean(axis=1), marker=markers[idx],
+                        plt.semilogy(x_idx, SeqValues, marker=markers[idx],
                                      ls='--', lw=2, color=colors[idx], label=name.split("_rep",1)[0])
                         
                         
-                        
-                        
-                plt.xticks(x_idx)
+                # Set the major and minor locators
+                ax.xaxis.set_major_locator(ticker.AutoLocator())
+                ax.xaxis.set_minor_locator(ticker.AutoMinorLocator())
+                ax.xaxis.grid(True, which='major', linestyle='-')
+                ax.xaxis.grid(True, which='minor', linestyle='--')
+                
+                ax.tick_params(axis='both', which='major', direction='in', width=3, length=10)
+                ax.tick_params(axis='both', which='minor', direction='in', width=2, length=8)
+                
                 plt.xlabel('Number of runs')
                 plt.ylabel(plotLabel)
                 plt.legend()