new figure for priors and improved labels for simulation results

simulation_scripts/plot_Q2_params.py

@@ -9,12 +9,13 @@
 from plotting import prepare_figure, draw_notches, draw_dot, globfigparams, colsgrey
 import numpy as np
 #pathoutrep = os.path.join(path, scenario, str(n), str(rep))
-horlims = {'sigmap':(0,0.0045),'mu':(0,0.022),'lambda':(0,0.09),'kappa':(0,0.4)}
-xticks = {'sigmap': (0, 0.003), 'mu':(0,0.02), 'lambda': (0,0.05), 'kappa': (0,0.2)}
+horlims = {'sigmap':(0,0.0045),'mu':(0,0.022),'lambda':(0,0.09),'kappa':(0,0.45)}
+#xticks = {'sigmap': (0, 0.003), 'mu':(0,0.02), 'lambda': (0,0.05), 'kappa': (0,0.2)}
+xticks = {'sigmap': (0, 0.002, 0.004), 'mu':(0,0.01,0.02), 'lambda': (0,0.04,0.08), 'kappa': (0,0.2,0.4)}
 skipbars = [('Q2kappa_base','kappa'),('Q2kappa_base','lambda'),('Q2kappa_base','mu'),('Q2kappa_lambdamu','kappa')]
 clipfirstvalues = ['mu','lambda','kappa']
 xlabels = {'sigmap': '$\\sigma$ [$\\mathrm{m}^3$ $\\mathrm{m}^{-3}$]', 'mu': '$\\mu$ [$\\mathrm{m}^3$ $\\mathrm{m}^{-3}$]', 'lambda': '$\\lambda$ [-]', 'kappa':'$\\kappa$ [-]'}
-coltitles = {'sigmap': 'Error $\\sigma$', 'lambda': 'Multiplicative $\\lambda$', 'kappa': 'Noise coeff. $\\kappa$', 'mu': 'Additive $\\mu$'}
+coltitles = {'sigmap': 'Error $\\sigma$', 'lambda': 'Sensitivity $\\lambda$', 'kappa': 'Noise coeff. $\\kappa$', 'mu': 'Additive $\\mu$'}
 rowlabels = {'Q2kappa_base': 'No $\\mu$, $\\lambda$, $\\kappa$', 'Q1kappa': 'Full model', 'Q2kappa_lambdamu': 'No $\\kappa$'}
 lw=0.5
 ms=5
@@ -93,11 +94,12 @@ def plot_Q2_params(parameters, scenarios = ['Q2kappa_base', 'Q2lambdamu_base', '
     axleg.barh(1, barl, height, color='#dddddd')
     draw_notches(axleg,biasl,1,height=height)
     draw_dot(axleg, dotl, 1)
-    axleg.annotate('Bias', xy=(biasl,1-0.5*height), xytext=(biasl-0.02,0.3),ha='center',va='top')
+    axleg.annotate('$b$', xy=(biasl,1-0.5*height), xytext=(biasl-0.02,0.3),ha='center',va='top')
     axleg.plot((biasl, biasl-0.01),(1-0.5*height-0.1, 0.35),lw=plt.rcParams['axes.linewidth'],c=globfigparams['fontcolour'])
     axleg.annotate('RMSE', xy=(barl,1+0.5*height), xytext=(barl+0.00,1.7),ha='center',va='bottom')
     axleg.plot((barl, barl),(1+0.5*height+0.1, 1.67),lw=plt.rcParams['axes.linewidth'],c=globfigparams['fontcolour'])
-    axleg.annotate('Posterior', xy=(dotl, 1-0.5*height), xytext=(dotl+0.1,0.3), ha='center', va='top')
+    axleg.annotate('$s_p$', xy=(dotl, 1-0.5*height), xytext=(dotl+0.1,0.3), ha='center', va='top')
+
     axleg.plot((dotl+0.02, dotl+0.04),(1-0.3, 0.35),lw=plt.rcParams['axes.linewidth'],c=globfigparams['fontcolour'])
 
     '''
@@ -117,7 +119,7 @@ def plot_Q2_params(parameters, scenarios = ['Q2kappa_base', 'Q2lambdamu_base', '
 
 def plot_Q2_params_new(parameters, scenarios = ['Q2kappa_base', 'Q2lambdamu_base', 'Q2kappa_lambdamu'], ns = [100,250,500], horlims = {}):
     ncols = len(parameters)    
-    fig, axs = prepare_figure(nrows=2, ncols=ncols, figsize_columns=(1.7, 0.75), sharex='col', sharey=True, bottom=0.16, left=0.11, right=0.85, top=0.82, hspace=0.5, wspace=0.25)
+    fig, axs = prepare_figure(nrows=2, ncols=ncols, figsize_columns=(1.7, 0.78), sharex=False, sharey=True, bottom=0.16, left=0.11, right=0.85, top=0.82, hspace=0.7, wspace=0.25)
 
     import input_output
     import os
@@ -152,7 +154,8 @@ def plot_Q2_params_new(parameters, scenarios = ['Q2kappa_base', 'Q2lambdamu_base
             axs[0,jparam].plot(horline,verticalpos,color=colsgrey[jn],lw=lw,alpha=0.75)
 
         if parameter in xticks:
-            axs[0,jparam].set_xticks(xticks[parameter])               
+            axs[0,jparam].set_xticks(xticks[parameter])              
+
         title=axs[0,jparam].set_title(coltitles[parameter],size=8)
         title.set_position((0.5,1.21))
         #axs[jpanel].set_title(titles[jpanel])
@@ -172,19 +175,17 @@ def plot_Q2_params_new(parameters, scenarios = ['Q2kappa_base', 'Q2lambdamu_base
         ydot = np.array([metrics[scenario][parameter][n][metricdot] for scenario in scenarios])
         #y1 = np.array([metrics[scenario][n][rowmetrics[1]] for n in ns])
         for jscenario,scenario in enumerate(scenarios):
-            if (scenario,parameter) not in skipbars:
-                axs[1,jparam].barh(verticalpos[jscenario], ybar[jscenario], height, color=colsgrey[2])
-            else:
-                pass
-                #axs[1,jparam].text(0.06,0,'$\\rightarrow$ 0',va='bottom',ha='left',transform=axs[1,jparam].transAxes)        
+            #if (scenario,parameter) not in skipbars:
+            axs[1,jparam].barh(verticalpos[jscenario], ybar[jscenario], height, color=colsgrey[2])
+
             draw_notches(axs[1,jparam],yline[jscenario],verticalpos[jscenario],height=height)
             draw_dot(axs[1,jparam],ydot[jscenario],verticalpos[jscenario])
         #axs[1,jparam].set_yticks([])
         if parameter in xticks:
             axs[1,jparam].set_xticks(xticks[parameter])
-        if parameter in xlabels:
+        if parameter in xlabels:            
             axs[1,jparam].set_xlabel(xlabels[parameter])
-
+        
     for jrow in np.arange(2):
         axs[jrow,0].set_yticks(verticalpos)
         axs[jrow,0].set_yticklabels([rowlabels[scenario] for scenario in scenarios])
@@ -205,23 +206,23 @@ def plot_Q2_params_new(parameters, scenarios = ['Q2kappa_base', 'Q2lambdamu_base
     axleg.set_yticks([])
     axleg.set_xlim((0,1))
     barl=0.25
-    biasl=0.07
+    biasl=0.09
     dotl=0.35
     axleg.barh(1, barl, height, color='#dddddd')
     draw_notches(axleg,biasl,1,height=height)
     draw_dot(axleg, dotl, 1)
-    axleg.annotate('Bias', xy=(biasl,1-0.5*height), xytext=(biasl-0.1,0.3),ha='center',va='top', size=fontsizelegend)
-    axleg.plot((biasl-0.01, biasl-0.025),(1-0.5*height-0.1, 0.35),lw=plt.rcParams['axes.linewidth'],c=globfigparams['fontcolour'])
+    axleg.annotate('$b$', xy=(biasl,1-0.5*height), xytext=(biasl-0.018,0.3),ha='center',va='top', size=fontsizelegend)
+    axleg.plot((biasl, biasl),(1-0.5*height-0.1, 0.35),lw=plt.rcParams['axes.linewidth'],c=globfigparams['fontcolour'])
     axleg.annotate('RMSE', xy=(barl,1+0.5*height), xytext=(barl+0.00,1.7),ha='center',va='bottom', size=fontsizelegend)
     axleg.plot((barl, barl),(1+0.5*height+0.1, 1.67),lw=plt.rcParams['axes.linewidth'],c=globfigparams['fontcolour'])
-    axleg.annotate('Posterior', xy=(dotl, 1-0.4*height), xytext=(dotl,0.3), ha='center', va='top', size=fontsizelegend)
+    axleg.annotate('$s_p$', xy=(dotl, 1-0.4*height), xytext=(dotl,0.3), ha='center', va='top', size=fontsizelegend)
     axleg.plot((dotl+0.00, dotl+0.00),(1-0.3, 0.35),lw=plt.rcParams['axes.linewidth'],c=globfigparams['fontcolour'])        
 
     yheader=1.07
     xheader=-0.7
     colhead = 'k'
     axs[0,0].text(xheader,yheader,'a) Dependence of RMSE on sample size',transform=axs[0,0].transAxes,va='bottom',ha='left', color=colhead)
-    axs[1,0].text(xheader,yheader,'b) Bias and posterior uncertainty for 500 samples',transform=axs[1,0].transAxes,va='bottom',ha='left', color=colhead)
+    axs[1,0].text(xheader,yheader,'b) Average bias and posterior standard deviation for 500 samples',transform=axs[1,0].transAxes,va='bottom',ha='left', color=colhead)
     axs[0,0].text(0.5,0.995,'\\textbf{Estimation accuracy in the simulation study: full and simplified models}', transform=fig.transFigure, color='k', va='top', ha='center')
     for ax in axs.flatten():
         ax.tick_params(axis='y', which='both',length=0)

simulation_scripts/plot_Q3_params.py

@@ -13,7 +13,7 @@
 clipfirstvalues = ['mu','lambda','kappa']
 scenariolabels = {'Q1kappa':'baseline','Q3dof':'prior tails', 'Q3priorfactor':'prior scale', 'Q3beta':'moisture marginal', 'Q3logisticspline':'moisture seasonal', 'Q3ar1':'autocorrelation', 'Q3studenttsim':'t error simulation', 'Q3studenttinference':'t error inference', 'Q3studenttsiminference':'t error'}
 xlabels = {'sigmap': '$\\sigma$ [$\\mathrm{m}^3$ $\\mathrm{m}^{-3}$]', 'mu': '$\\mu$ [$\\mathrm{m}^3$ $\\mathrm{m}^{-3}$]', 'lambda': '$\\lambda$ [-]', 'kappa':'$\\kappa$ [-]'}
-coltitles = {'sigmap': 'Error $\\sigma$', 'lambda': 'Multiplicative $\\lambda$', 'kappa': 'Noise coeff. $\\kappa$', 'mu': 'Additive $\\mu$'}
+coltitles = {'sigmap': 'Error $\\sigma$', 'lambda': 'Sensitivity $\\lambda$', 'kappa': 'Noise coeff. $\\kappa$', 'mu': 'Offset $\\mu$'}
 def plot_Q3_params(parameters, scenarios = ['Q3dof', 'Q3priorfactor', 'Q3beta', 'Q3logisticspline', 'Q3ar1', 'Q3studenttsim', 'Q3studenttinference', 'Q3studenttsiminference'], n=500, horlims = {}):
     ncols = len(parameters)
     sharex='col'

simulation_scripts/plot_priors.py

@@ -0,0 +1,45 @@
+'''
+Created on Feb 27, 2018
+
+@author: zwieback
+'''
+# sigma: exponential: mean 0.1^2
+#m and mu: T(0,0.3^2;4)
+#l : T(1,0.3^2;4)
+import numpy as np
+from scipy.stats import t, norm
+import matplotlib.pyplot as plt
+from plotting import colsgrey,prepare_figure,globfigparams
+
+
+labels = ['$\\sigma^2$ $[(\\mathrm{m}^3 \\, \\mathrm{m}^{-3})^2]$', '$m$ and $\\mu$ [$\\mathrm{m}^3 \\, \\mathrm{m}^{-3}$]', '$l$ [-]', '$\\lambda$ [-]']
+
+fig, axs = prepare_figure(nrows=1, ncols=4, figsize_columns=(1.9, 0.5), sharex=False, sharey=False, bottom=0.22, left=0.02, right=0.985, top=0.92, hspace=0.39)
+lw=1.0
+col = 'k'
+xsigma2 = np.linspace(0,0.5)
+axs[0].plot(xsigma2,np.exp(-xsigma2/(0.1)), lw=lw, c=col)
+
+xm = np.linspace(-0.7,0.7)
+df=4
+scale = 0.3/t.std(df=df)
+dist = t(df=df,scale=scale)
+axs[1].plot(xm, dist.pdf(xm), lw=lw, c=col)
+axs[3].plot(xm, dist.pdf(xm), lw=lw, c=col)
+
+xl = np.linspace(0,2.0)
+df=4
+scale = 0.3/t.std(df=df)
+dist = t(df=df,scale=scale,loc=1.0)
+axs[2].plot(xl, dist.pdf(xl), lw=lw, c=col)
+
+for j,ax in enumerate(axs.flatten()):
+    for loc, spine in ax.spines.items():
+        if loc in ['right','top','left']:
+            spine.set_color('none')    
+            ax.set_yticks([])
+            ylim=(0,ax.get_ylim()[1])
+            ax.set_ylim(ylim)           
+            ax.text(0.5,-0.2,labels[j],ha='center', va='top', transform=ax.transAxes)#, color=globfigparams['fontcolour'])
+ax.text(0.5,0.98,'Prior distributions',ha='center',va='top', transform=fig.transFigure, size=8, color='k')            
+plt.savefig('C:\\Work\\SMAP\\simulations\\priors.pdf')