diff --git a/graphics.py b/graphics.py
index 8c5dca07103f40a4b88e7ae8079e45cbd4c97877..64e5227137b2d25fe6de532916bd70a5e7f294dd 100644
--- a/graphics.py
+++ b/graphics.py
@@ -133,6 +133,7 @@ def plot_CBL_PE(exp,figtitle,parameter_id=0,plot_spread=False,ax=None):
         parameter_number = exp.nr.parameter_number
         par_tran = exp.da.backgrounds[:,-parameter_number:,:][:,parameter_id,:]
         par_phys = pt[parameter_id](par_tran,kind='inv')
+        initpar_phys = pt[parameter_id](exp.da.initial_perturbed_parameters,kind='inv')
         t = exp.da.time/3600.
 
         plot_parameter_histogram = True
@@ -187,7 +188,18 @@ def plot_CBL_PE(exp,figtitle,parameter_id=0,plot_spread=False,ax=None):
         else:
             ax1 = ax
 
-        ax1.step(t,par_phys.mean(axis=1),color=mean_color)
+        # Initial parameter distribution
+        ax1.step([0,t[0]],np.median(initpar_phys,axis=1)*np.array([1,1]),color=mean_color)
+        bmin = np.percentile(initpar_phys, 10)
+        bmax = np.percentile(initpar_phys, 90)
+        ax1.fill_between([0,t[0]],\
+                    y1=[bmin,bmin],\
+                    y2=[bmax,bmax],\
+                    color=spread_color,
+                    edgecolor=None,
+                    alpha=0.3)
+        # Later times
+        ax1.step(t,np.median(par_phys,axis=1),color=mean_color)
         for i in range(t.size-1):
             bmin = np.percentile(par_phys, 10, axis=1)[i+1]
             bmax = np.percentile(par_phys, 90, axis=1)[i+1]