From bd905360e0be1ea8584452d5e84637d0877d3bd1 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 8 Sep 2023 10:26:21 +0200 Subject: [PATCH] Slice: provide convergence diagnostics even for low number of draws --- matlab/convergence_diagnostics/mcmc_diagnostics.m | 15 ++++++++++----- .../mcmc_diagnostics_core.m | 4 ++-- 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/matlab/convergence_diagnostics/mcmc_diagnostics.m b/matlab/convergence_diagnostics/mcmc_diagnostics.m index 8ad512a68..1a3bc0cbb 100644 --- a/matlab/convergence_diagnostics/mcmc_diagnostics.m +++ b/matlab/convergence_diagnostics/mcmc_diagnostics.m @@ -110,7 +110,7 @@ update_last_mh_history_file(MetropolisFolder, ModelName, record); PastDraws = sum(record.MhDraws,1); NumberOfDraws = PastDraws(1); -if NumberOfDraws<=2000 +if ~strcmp(options_.posterior_sampler_options.posterior_sampling_method,'slice') && NumberOfDraws<=2000 warning('MCMC convergence diagnostics are not computed because the total number of iterations is not bigger than 2000!'); return end @@ -210,8 +210,13 @@ if nblck == 1 % Brooks and Gelman tests need more than one block return end -Origin = 1000; -StepSize = ceil((NumberOfDraws-Origin)/100);% So that the computational time does not +if strcmp(options_.posterior_sampler_options.posterior_sampling_method,'slice') && NumberOfDraws<2000 + Origin = 1; + StepSize = 1; +else + Origin = 1000; + StepSize = ceil((NumberOfDraws-Origin)/100);% So that the computational time does not +end ALPHA = 0.2; % increase too much with the number of simulations. time = 1:NumberOfDraws; xx = Origin:StepSize:NumberOfDraws; @@ -421,9 +426,9 @@ for iter = Origin:StepSize:NumberOfDraws ligne = ligne+1; linea = ceil(options_.mh_drop*iter); n = iter-linea+1; - cinf = round(n*ALPHA/2); + cinf = max(1,round(n*ALPHA/2)); csup = round(n*(1-ALPHA/2)); - CINF = round(nblck*n*ALPHA/2); + CINF = max(1,round(nblck*n*ALPHA/2)); CSUP = round(nblck*n*(1-ALPHA/2)); temp = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2); MDIAG(ligne,1) = temp(CSUP,1)-temp(CINF,1); diff --git a/matlab/convergence_diagnostics/mcmc_diagnostics_core.m b/matlab/convergence_diagnostics/mcmc_diagnostics_core.m index a9ab4581e..fb0dfd32c 100644 --- a/matlab/convergence_diagnostics/mcmc_diagnostics_core.m +++ b/matlab/convergence_diagnostics/mcmc_diagnostics_core.m @@ -112,9 +112,9 @@ for j=fpar:npar window_iter = window_iter+1; linea = ceil(mh_drop*iter); %compute first non-discarded draw; drops fraction of sample at each iteration for computational efficiency, see Brooks/Gelman (1998), p.438 n = iter-linea+1; %number of draws from each block in current batch - cinf = round(n*ALPHA/2); %lower bound for alpha percentile of within series + cinf = max(1,round(n*ALPHA/2)); %lower bound for alpha percentile of within series csup = round(n*(1-ALPHA/2)); %upper bound for alpha percentile of within series - CINF = round(nblck*n*ALPHA/2); %lower bound for alpha percentile of pooled series + CINF = max(1,round(nblck*n*ALPHA/2)); %lower bound for alpha percentile of pooled series CSUP = round(nblck*n*(1-ALPHA/2)); %upper bound for alpha percentile of pooled series temp = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2); %extract pooled draws in current batch UDIAG(window_iter,1,j-fpar+1) = temp(CSUP,1)-temp(CINF,1); %length of total sequence interval