added code to plot parameter distributions

graphics.py

@@ -131,10 +131,47 @@ def plot_CBL_PE(exp,figtitle,parameter_id=0,plot_spread=False,ax=None):
 
         pt = exp.nr.parameter_transform
         parameter_number = exp.nr.parameter_number
-        b = exp.da.backgrounds[:,-parameter_number:,:][:,parameter_id,:]
-        par_b = pt[parameter_id](exp.da.backgrounds[:,-parameter_number:,:][:,parameter_id,:],kind='inv')
+        par_tran = exp.da.backgrounds[:,-parameter_number:,:][:,parameter_id,:]
+        par_phys = pt[parameter_id](par_tran,kind='inv')
         t = exp.da.time/3600.
 
+        plot_parameter_histogram = True
+        if plot_parameter_histogram:
+            print(par_tran[-1,:])
+            p10beg_tran = np.percentile(par_tran[0,:],  10)
+            p90beg_tran = np.percentile(par_tran[0,:],  90)
+            p10end_tran = np.percentile(par_tran[-1,:], 10)
+            p90end_tran = np.percentile(par_tran[-1,:], 90)
+            p10beg_phys = np.percentile(par_phys[0,:],  10)
+            p90beg_phys = np.percentile(par_phys[0,:],  90)
+            p10end_phys = np.percentile(par_phys[-1,:], 10)
+            p90end_phys = np.percentile(par_phys[-1,:], 90)
+            #
+            fig = p.figure(151)
+            fig.set_size_inches(6,3)
+            #
+            ax1 = fig.add_subplot(1,2,1)
+            ax1.hist(par_phys[0,:], bins=np.linspace(0,5,51), density=True, histtype='step', color='blue', label='initial')
+            ax1.axvspan(p10beg_phys,p90beg_phys,color='blue',lw=0,alpha=0.1)
+            ax1.hist(par_phys[-1,:], bins=np.linspace(0,5,51), density=True, histtype='step', color='red', label='final')
+            ax1.axvspan(p10end_phys,p90end_phys,color='red',lw=0,alpha=0.1)
+            ax1.axvline(x=exp.nr.parameter_true,color='black',dashes=[3,1])
+            ax1.set_xlabel(r'$p$')
+            ax1.set_ylabel('probability density')
+            #
+            ax2 = fig.add_subplot(1,2,2)
+            ax2.hist(par_tran[0,:], bins=np.linspace(np.log(0.01),np.log(5),51), density=True, histtype='step', color='blue', label='initial')
+            ax2.axvspan(p10beg_tran,p90beg_tran,color='blue',lw=0,alpha=0.1)
+            ax2.hist(par_tran[-1,:], bins=np.linspace(np.log(0.01),np.log(5),51), density=True, histtype='step', color='red', label='final')
+            ax2.axvspan(p10end_tran,p90end_tran,color='red',lw=0,alpha=0.1)
+            ax2.axvline(x=pt[parameter_id](exp.nr.parameter_true,kind='dir'),color='black',dashes=[3,1])
+            ax2.set_xlabel(r'log($p$)')
+            ax2.set_ylabel('probability density')
+            #
+            fig.tight_layout()
+            fig.savefig('p_histogram_%s.png'%exp.label,format='png',dpi=300)
+            p.close(fig)
+
         if ax is None:
             printfigure = True
         else:
@@ -151,10 +188,10 @@ def plot_CBL_PE(exp,figtitle,parameter_id=0,plot_spread=False,ax=None):
         else:
             ax1 = ax
 
-        ax1.step(t,par_b.mean(axis=1),color=mean_color)
+        ax1.step(t,par_phys.mean(axis=1),color=mean_color)
         for i in range(t.size-1):
-            bmin = np.percentile(par_b, 10, axis=1)[i+1]
-            bmax = np.percentile(par_b, 90, axis=1)[i+1]
+            bmin = np.percentile(par_phys, 10, axis=1)[i+1]
+            bmax = np.percentile(par_phys, 90, axis=1)[i+1]
             ax1.fill_between(t[i:i+2],\
                         y1=[bmin,bmin],\
                         y2=[bmax,bmax],\
@@ -166,8 +203,8 @@ def plot_CBL_PE(exp,figtitle,parameter_id=0,plot_spread=False,ax=None):
 
         if printfigure and plot_spread:
             ax2 = fig.add_subplot(ii,1,2)
-            ax2.plot(t,par_b.std(axis=1),color=mean_color,label='physical')
-            ax2.plot(t,b.std(axis=1),color=mean_color_2,label='transformed', linestyle='--')
+            ax2.plot(t,par_phys.std(axis=1),color=mean_color,label='physical')
+            ax2.plot(t,par_tran.std(axis=1),color=mean_color_2,label='transformed', linestyle='--')
             ax2.legend()
             ax2.set_title('Parameter spread')