From bf674fa59fc300f6a5ba3aaf30eecfbff2521230 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx,de>
Date: Thu, 6 Mar 2014 08:10:51 +0100
Subject: [PATCH] Deal with cases where too few draws are available for
 posterior moments

---
 matlab/GetPosteriorParametersStatistics.m |  8 +++++
 matlab/posterior_moments.m                | 43 +++++++++++++----------
 2 files changed, 33 insertions(+), 18 deletions(-)

diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index f92c2629e..cbcc318b4 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/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
diff --git a/matlab/posterior_moments.m b/matlab/posterior_moments.m
index ff9eeb8e1..7e122cf46 100644
--- a/matlab/posterior_moments.m
+++ b/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;
+
+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
-[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)]);
-
 density = [];
 if info
     number_of_grid_points = 2^9;      % 2^9 = 512 !... Must be a power of two.