Deal with cases where too few draws are available for posterior moments

(cherry picked from commit bf674fa5)

matlab/GetPosteriorParametersStatistics.m

@@ -75,6 +75,14 @@ catch
     [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);
     disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
 end
+num_draws=length(NumberOfDraws*options_.mh_nblck);
+hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
+if hpd_draws<2
+    fprintf('posterior_moments: There are not enough draws computes to compute HPD Intervals. Skipping their computation.\n')
+end
+if num_draws<9
+    fprintf('posterior_moments: There are not enough draws computes to compute deciles. Skipping their computation.\n')
+end
 if np
     type = 'parameters';
     if TeX

matlab/posterior_moments.m

@@ -48,25 +48,32 @@ post_var = var(xx);
 
 number_of_draws = length(xx);
 hpd_draws = round((1-mh_conf_sig)*number_of_draws);
-kk = zeros(hpd_draws,1);
-jj = number_of_draws-hpd_draws;
-for ii = 1:hpd_draws
-    kk(ii) = xx(jj)-xx(ii);
-    jj = jj + 1;
-end
-[kmin,idx] = min(kk);
-hpd_interval = [xx(idx) xx(idx)+kmin];
-
-post_deciles = xx([round(0.1*number_of_draws) ...
-                   round(0.2*number_of_draws) ...
-                   round(0.3*number_of_draws) ...
-                   round(0.4*number_of_draws) ...
-                   round(0.5*number_of_draws) ...
-                   round(0.6*number_of_draws) ...
-                   round(0.7*number_of_draws) ...
-                   round(0.8*number_of_draws) ...
-                   round(0.9*number_of_draws)]);
 
+if hpd_draws>2
+    kk = zeros(hpd_draws,1);
+    jj = number_of_draws-hpd_draws;
+    for ii = 1:hpd_draws
+        kk(ii) = xx(jj)-xx(ii);
+        jj = jj + 1;
+    end
+    [kmin,idx] = min(kk);
+    hpd_interval = [xx(idx) xx(idx)+kmin];
+else
+    hpd_interval=NaN(1,2);    
+end
+if length(xx)>9
+    post_deciles = xx([round(0.1*number_of_draws) ...
+                       round(0.2*number_of_draws) ...
+                       round(0.3*number_of_draws) ...
+                       round(0.4*number_of_draws) ...
+                       round(0.5*number_of_draws) ...
+                       round(0.6*number_of_draws) ...
+                       round(0.7*number_of_draws) ...
+                       round(0.8*number_of_draws) ...
+                       round(0.9*number_of_draws)]);
+else
+    post_deciles=NaN(9,1);
+end
 density = [];
 if info
     number_of_grid_points = 2^9;      % 2^9 = 512 !... Must be a power of two.