diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 4434a2666c7b247d083386cda61e9e565aa643dd..c4512438316d5e3ba34b387db5d4d179ab7cddd8 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -714,6 +714,9 @@ options_.convergence.geweke.geweke_interval=[0.2 0.5];
 options_.convergence.rafterylewis.indicator=false;
 options_.convergence.rafterylewis.qrs=[0.025 0.005 0.95];
 
+%tolerance for Modified Harmonic Mean estimator
+options_.marginal_data_density.harmonic_mean.tolerance = 0.01;
+
 % Options for lmmcp solver
 options_.lmmcp.status = false;
 
diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m
index 1e87281610b6dd3a35ca15ddf34f70e338929dac..1179fb50c1dd6f7062d1af2649b28341fdc25e1e 100644
--- a/matlab/marginal_density.m
+++ b/matlab/marginal_density.m
@@ -100,7 +100,7 @@ while check_coverage
         marginal(linee,:) = [p, lpost_mode-log(tmp/((TotalNumberOfMhDraws-TODROP)*nblck))];
         warning(warning_old_state);
     end
-    if abs((marginal(9,2)-marginal(1,2))/marginal(9,2)) > 0.01 || isinf(marginal(1,2))
+    if abs((marginal(9,2)-marginal(1,2))/marginal(9,2)) > options_.marginal_data_density.harmonic_mean.tolerance || isinf(marginal(1,2))
         fprintf('\n')
         if increase == 1
             disp('Estimation::marginal density: The support of the weighting density function is not large enough...')