stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

fixed for the case when mcmc is incomplete WITHIN a block file (useful for expensive models and expensive methods like slice or TaRB)

mr#2177
Marco Ratto 2023-11-10 17:51:41 +01:00
parent 53b57da8ba
commit f102a992aa
2 changed files with 34 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
matlab/posterior_sampler_core.m
View File

@ -141,31 +141,32 @@ for curr_block = fblck:nblck
options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block);
end
mh_recover_flag=0;
if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2 && OpenOldFile(curr_block)
% this should be done whatever value of load_mh_file
load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
draw_iter = size(neval_this_chain,2)+1;
draw_index_current_file = draw_iter+fline(curr_block)-1;
feval_this_chain = sum(sum(neval_this_chain));
feval_this_file = sum(sum(neval_this_chain));
if feval_this_chain>draw_index_current_file-fline(curr_block)
% non Metropolis type of sampler
accepted_draws_this_chain = draw_index_current_file-fline(curr_block);
accepted_draws_this_file = draw_index_current_file-fline(curr_block);
else
accepted_draws_this_chain = 0;
accepted_draws_this_file = 0;
end
mh_recover_flag=1;
set_dynare_random_generator_state(LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal);
last_draw(curr_block,:)=x2(draw_index_current_file-1,:);
last_posterior(curr_block)=logpo2(draw_index_current_file-1);
OpenOldFile(curr_block) = 0;
else
if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2
load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
draw_iter = size(neval_this_chain,2)+1;
draw_index_current_file = draw_iter;
feval_this_chain = sum(sum(neval_this_chain));
feval_this_file = sum(sum(neval_this_chain));
if feval_this_chain>draw_iter-1
% non Metropolis type of sampler
accepted_draws_this_chain = draw_iter-1;
accepted_draws_this_file = draw_iter-1;
else
accepted_draws_this_chain = 0;
accepted_draws_this_file = 0;
end
mh_recover_flag=1;
set_dynare_random_generator_state(LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal);
last_draw(curr_block,:)=x2(draw_iter-1,:);
last_posterior(curr_block)=logpo2(draw_iter-1);
if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
OpenOldFile(curr_block) = 0;
else
x2 = zeros(InitSizeArray(curr_block),npar);

12
matlab/posterior_sampler_initialization.m
View File

@ -465,7 +465,7 @@ elseif options_.mh_recover
AllMhFiles = dir([BaseName '_mh*_blck*.mat']);
TotalNumberOfMhFiles = length(AllMhFiles)-length(dir([BaseName '_mh_tmp*_blck*.mat']));
% Quit if no crashed mcmc chain can be found as there are as many files as expected
if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
if (ExpectedNumberOfMhFilesPerBlock>LastFileNumberInThePreviousMh) && (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
if isnumeric(options_.parallel)
fprintf('%s: It appears that you don''t need to use the mh_recover option!\n',dispString);
fprintf(' You have to edit the mod file and remove the mh_recover option\n');
@ -476,15 +476,23 @@ elseif options_.mh_recover
% 2. Something needs to be done; find out what
% Count the number of saved mh files per block.
NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1);
is_chain_complete = true(NumberOfBlocks,1);
for b = 1:NumberOfBlocks
BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']);
NumberOfMhFilesPerBlock(b) = length(BlckMhFiles)-length(dir([BaseName '_mh_tmp*_blck' int2str(b) '.mat']));
if (ExpectedNumberOfMhFilesPerBlock==LastFileNumberInThePreviousMh)
% here we need to check for the number of lines
tmpdata = load([BaseName '_mh' int2str(LastFileNumberInThePreviousMh) '_blck' int2str(b) '.mat'],'logpo2');
if length(tmpdata.logpo2) == LastLineNumberInThePreviousMh
is_chain_complete(b) = false;
end
end
end
% Find FirstBlock (First block), an integer targeting the crashed mcmc chain.
FirstBlock = 1; %initialize
FBlock = zeros(NumberOfBlocks,1);
while FirstBlock <= NumberOfBlocks
if NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock
if (NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock) || not(is_chain_complete(FirstBlock))
fprintf('%s: Chain %u is not complete!\n', dispString,FirstBlock);
FBlock(FirstBlock)=1;
else