diff --git a/matlab/metropolis_hastings_initialization.m b/matlab/metropolis_hastings_initialization.m deleted file mode 100644 index d3d9bc93f..000000000 --- a/matlab/metropolis_hastings_initialization.m +++ /dev/null @@ -1,440 +0,0 @@ -function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d ] = ... - metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) -% function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d ] = ... -% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) -% Metropolis-Hastings initialization. -% -% INPUTS -% o TargetFun [char] string specifying the name of the objective -% function (posterior kernel). -% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values). -% o vv [double] (p*p) matrix, posterior covariance matrix (at the mode). -% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters. -% o dataset_ data structure -% o dataset_info dataset info structure -% o options_ options structure -% o M_ model structure -% o estim_params_ estimated parameters structure -% o bayestopt_ estimation options structure -% o oo_ outputs structure -% -% OUTPUTS -% o ix2 [double] (NumberOfBlocks*npar) vector of starting points for different chains -% o ilogpo2 [double] (NumberOfBlocks*1) vector of initial posterior values for different chains -% o ModelName [string] name of the mod-file -% o MetropolisFolder [string] path to the Metropolis subfolder -% o FirstBlock [scalar] number of the first MH chain to be run (not equal to 1 in case of recovery) -% o FirstLine [double] (NumberOfBlocks*1) vector of first draw in each chain (not equal to 1 in case of recovery) -% o npar [scalar] number of parameters estimated -% o NumberOfBlocks [scalar] Number of MCM chains requested -% o nruns [double] (NumberOfBlocks*1) number of draws in each chain -% o NewFile [scalar] (NumberOfBlocks*1) vector storing the number -% of the first MH-file to created for each chain when saving additional -% draws -% o MAX_nruns [scalar] maximum number of draws in each MH-file on the harddisk -% o d [double] (p*p) matrix, Cholesky decomposition of the posterior covariance matrix (at the mode). -% -% SPECIAL REQUIREMENTS -% None. - -% Copyright © 2006-2017 Dynare Team -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see . - -%Initialize outputs -ix2 = []; -ilogpo2 = []; -ModelName = []; -MetropolisFolder = []; -FirstBlock = []; -FirstLine = []; -npar = []; -NumberOfBlocks = []; -nruns = []; -NewFile = []; -MAX_nruns = []; -d = []; - -ModelName = M_.fname; -if ~isempty(M_.bvar) - ModelName = [ModelName '_bvar']; -end - -MetropolisFolder = CheckPath('metropolis',M_.dname); -BaseName = [MetropolisFolder filesep ModelName]; - -NumberOfBlocks = options_.mh_nblck; -nruns = ones(NumberOfBlocks,1)*options_.mh_replic; -npar = length(xparam1); -MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); -d = chol(vv); - -if ~options_.load_mh_file && ~options_.mh_recover - % Here we start a new Metropolis-Hastings, previous draws are discarded. - if NumberOfBlocks > 1 - disp('Estimation::mcmc: Multiple chains mode.') - else - disp('Estimation::mcmc: One Chain mode.') - end - % Delete old mh files if any... - files = dir([BaseName '_mh*_blck*.mat']); - if length(files) - delete([BaseName '_mh*_blck*.mat']); - disp('Estimation::mcmc: Old mh-files successfully erased!') - end - % Delete old Metropolis log file. - file = dir([ MetropolisFolder '/metropolis.log']); - if length(file) - delete([ MetropolisFolder '/metropolis.log']); - disp('Estimation::mcmc: Old metropolis.log file successfully erased!') - disp('Estimation::mcmc: Creation of a new metropolis.log file.') - end - fidlog = fopen([MetropolisFolder '/metropolis.log'],'w'); - fprintf(fidlog,'%% MH log file (Dynare).\n'); - fprintf(fidlog,['%% ' datestr(now,0) '.\n']); - fprintf(fidlog,' \n\n'); - fprintf(fidlog,'%% Session 1.\n'); - fprintf(fidlog,' \n'); - fprintf(fidlog,[' Number of blocks...............: ' int2str(NumberOfBlocks) '\n']); - fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']); - fprintf(fidlog,' \n'); - % Find initial values for the NumberOfBlocks chains... - if NumberOfBlocks > 1% Case 1: multiple chains - set_dynare_seed('default'); - fprintf(fidlog,[' Initial values of the parameters:\n']); - disp('Estimation::mcmc: Searching for initial values...') - ix2 = zeros(NumberOfBlocks,npar); - ilogpo2 = zeros(NumberOfBlocks,1); - for j=1:NumberOfBlocks - validate = 0; - init_iter = 0; - trial = 1; - while validate == 0 && trial <= 10 - candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar); - if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub) - ix2(j,:) = candidate; - ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_); - if ~isfinite(ilogpo2(j)) % if returned log-density is - % Inf or Nan (penalized value) - validate = 0; - else - fprintf(fidlog,[' Blck ' int2str(j) ':\n']); - for i=1:length(ix2(1,:)) - fprintf(fidlog,[' params:' int2str(i) ': ' num2str(ix2(j,i)) '\n']); - end - fprintf(fidlog,[' logpo2: ' num2str(ilogpo2(j)) '\n']); - j = j+1; - validate = 1; - end - end - init_iter = init_iter + 1; - if init_iter > 100 && validate == 0 - disp(['Estimation::mcmc: I couldn''t get a valid initial value in 100 trials.']) - if options_.nointeractive - disp(['Estimation::mcmc: I reduce mh_init_scale by 10 percent:']) - options_.mh_init_scale = .9*options_.mh_init_scale; - disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.',options_.mh_init_scale)) - else - disp(['Estimation::mcmc: You should reduce mh_init_scale...']) - disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is equal to %f.',options_.mh_init_scale)) - options_.mh_init_scale = input('Estimation::mcmc: Enter a new value... '); - end - trial = trial+1; - end - end - if trial > 10 && ~validate - disp(['Estimation::mcmc: I''m unable to find a starting value for block ' int2str(j)]) - fclose(fidlog); - return - end - end - fprintf(fidlog,' \n'); - disp('Estimation::mcmc: Initial values found!') - skipline() - else% Case 2: one chain (we start from the posterior mode) - fprintf(fidlog,[' Initial values of the parameters:\n']); - candidate = transpose(xparam1(:));% - if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub) - ix2 = candidate; - ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_); - disp('Estimation::mcmc: Initialization at the posterior mode.') - skipline() - fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']); - for i=1:length(ix2(1,:)) - fprintf(fidlog,[' ' int2str(i) ':' num2str(ix2(1,i)) '\n']); - end - fprintf(fidlog,[' Blck ' int2str(1) 'logpo2:' num2str(ilogpo2) '\n']); - else - disp('Estimation::mcmc: Initialization failed...') - disp('Estimation::mcmc: The posterior mode lies outside the prior bounds.') - fclose(fidlog); - return - end - fprintf(fidlog,' \n'); - end - fprintf(fidlog,' \n'); - FirstBlock = 1; - FirstLine = ones(NumberOfBlocks,1); - NewFile = ones(NumberOfBlocks,1); - % Delete the mh-history files - delete_mh_history_files(MetropolisFolder, ModelName); - % Create a new record structure - fprintf(['Estimation::mcmc: Write details about the MCMC... ']); - AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns); - AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns; - record.Nblck = NumberOfBlocks; - record.MhDraws = zeros(1,3); - record.MhDraws(1,1) = nruns(1); - record.MhDraws(1,2) = AnticipatedNumberOfFiles; - record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile; - record.MAX_nruns=MAX_nruns; - record.AcceptanceRatio = zeros(1,NumberOfBlocks); - for j=1:NumberOfBlocks - % we set a different seed for the random generator for each block then we record the corresponding random generator state (vector) - set_dynare_seed(options_.DynareRandomStreams.seed+j); - % record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name - [record.InitialSeeds(j).Unifor,record.InitialSeeds(j).Normal] = get_dynare_random_generator_state(); - end - record.InitialParameters = ix2; - record.InitialLogPost = ilogpo2; - record.LastParameters = zeros(NumberOfBlocks,npar); - record.LastLogPost = zeros(NumberOfBlocks,1); - record.LastFileNumber = AnticipatedNumberOfFiles ; - record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile; - record.MCMCConcludedSuccessfully = 0; - fprintf('Ok!\n'); - id = write_mh_history_file(MetropolisFolder, ModelName, record); - disp(['Estimation::mcmc: Details about the MCMC are available in ' BaseName '_mh_history_' num2str(id) '.mat']) - skipline() - fprintf(fidlog,[' CREATION OF THE MH HISTORY FILE!\n\n']); - fprintf(fidlog,[' Expected number of files per block.......: ' int2str(AnticipatedNumberOfFiles) '.\n']); - fprintf(fidlog,[' Expected number of lines in the last file: ' int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']); - fprintf(fidlog,['\n']); - for j = 1:NumberOfBlocks - fprintf(fidlog,[' Initial state of the Gaussian random number generator for chain number ',int2str(j),':\n']); - for i=1:length(record.InitialSeeds(j).Normal) - fprintf(fidlog,[' ' num2str(record.InitialSeeds(j).Normal(i)') '\n']); - end - fprintf(fidlog,[' Initial state of the Uniform random number generator for chain number ',int2str(j),':\n']); - for i=1:length(record.InitialSeeds(j).Unifor) - fprintf(fidlog,[' ' num2str(record.InitialSeeds(j).Unifor(i)') '\n']); - end - end - fprintf(fidlog,' \n'); - fclose(fidlog); -elseif options_.load_mh_file && ~options_.mh_recover - % Here we consider previous mh files (previous mh did not crash). - disp('Estimation::mcmc: I am loading past Metropolis-Hastings simulations...') - load_last_mh_history_file(MetropolisFolder, ModelName); - if ~isnan(record.MCMCConcludedSuccessfully) && ~record.MCMCConcludedSuccessfully - error('Estimation::mcmc: You are trying to load an MCMC that did not finish successfully. Please use mh_recover.') - end - record.MCMCConcludedSuccessfully=0; %reset indicator for this run - mh_files = dir([ MetropolisFolder filesep ModelName '_mh*.mat']); - if ~length(mh_files) - error('Estimation::mcmc: I cannot find any MH file to load here!') - end - fidlog = fopen([MetropolisFolder '/metropolis.log'],'a'); - fprintf(fidlog,'\n'); - fprintf(fidlog,['%% Session ' int2str(length(record.MhDraws(:,1))+1) '.\n']); - fprintf(fidlog,' \n'); - fprintf(fidlog,[' Number of blocks...............: ' int2str(NumberOfBlocks) '\n']); - fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']); - fprintf(fidlog,' \n'); - past_number_of_blocks = record.Nblck; - if past_number_of_blocks ~= NumberOfBlocks - disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!') - disp(['Estimation::mcmc: You declared ' int2str(NumberOfBlocks) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.']) - disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ]) - NumberOfBlocks = past_number_of_blocks; - options_.mh_nblck = NumberOfBlocks; - end - % I read the last line of the last mh-file for initialization of the new Metropolis-Hastings simulations: - LastFileNumber = record.LastFileNumber; - LastLineNumber = record.LastLineNumber; - if LastLineNumber < MAX_nruns - NewFile = ones(NumberOfBlocks,1)*LastFileNumber; - FirstLine = ones(NumberOfBlocks,1)*(LastLineNumber+1); - else - NewFile = ones(NumberOfBlocks,1)*(LastFileNumber+1); - FirstLine = ones(NumberOfBlocks,1); - end - ilogpo2 = record.LastLogPost; - ix2 = record.LastParameters; - FirstBlock = 1; - NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1); - fprintf('Estimation::mcmc: I am writing a new mh-history file... '); - record.MhDraws = [record.MhDraws;zeros(1,3)]; - NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber; - NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile; - AnticipatedNumberOfFiles = ceil(NumberOfDrawsToBeSaved/MAX_nruns); - AnticipatedNumberOfLinesInTheLastFile = NumberOfDrawsToBeSaved - (AnticipatedNumberOfFiles-1)*MAX_nruns; - record.LastFileNumber = LastFileNumber + AnticipatedNumberOfFiles; - record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile; - record.MhDraws(end,1) = nruns(1); - record.MhDraws(end,2) = AnticipatedNumberOfFiles; - record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile; - record.MAX_nruns = [record.MAX_nruns;MAX_nruns]; - record.InitialSeeds = record.LastSeeds; - write_mh_history_file(MetropolisFolder, ModelName, record); - fprintf('Done.\n') - disp(['Estimation::mcmc: Ok. I have loaded ' int2str(NumberOfPreviousSimulations) ' simulations.']) - skipline() - fclose(fidlog); -elseif options_.mh_recover - % The previous metropolis-hastings crashed before the end! I try to recover the saved draws... - disp('Estimation::mcmc: Recover mode!') - load_last_mh_history_file(MetropolisFolder, ModelName); - NumberOfBlocks = record.Nblck;% Number of "parallel" mcmc chains. - options_.mh_nblck = NumberOfBlocks; - - %% check consistency of options - if record.MhDraws(end,1)~=options_.mh_replic - fprintf('\nEstimation::mcmc: You cannot specify a different mh_replic than in the chain you are trying to recover\n') - fprintf('Estimation::mcmc: I am resetting mh_replic to %u\n',record.MhDraws(end,1)) - options_.mh_replic=record.MhDraws(end,1); - nruns = ones(NumberOfBlocks,1)*options_.mh_replic; - end - - if ~isnan(record.MAX_nruns(end,1)) %field exists - if record.MAX_nruns(end,1)~=MAX_nruns - fprintf('\nEstimation::mcmc: You cannot specify a different MaxNumberOfBytes than in the chain you are trying to recover\n') - fprintf('Estimation::mcmc: I am resetting MAX_nruns to %u\n',record.MAX_nruns(end,1)) - MAX_nruns=record.MAX_nruns(end,1); - end - end - %% do tentative initialization as if full last run of MCMC were to be redone - if size(record.MhDraws,1) == 1 - OldMhExists = 0;% The crashed metropolis was the first session. - else - OldMhExists = 1;% The crashed metropolis wasn't the first session. - end - % Default initialization: - if OldMhExists - ilogpo2 = record.LastLogPost; - ix2 = record.LastParameters; - else - ilogpo2 = record.InitialLogPost; - ix2 = record.InitialParameters; - end - % Set NewFile, a NumberOfBlocks*1 vector of integers, and FirstLine (first line), a NumberOfBlocks*1 vector of integers. - % Relevant for starting at the correct file and potentially filling a file from the previous run that was not yet full - if OldMhExists - LastLineNumberInThePreviousMh = record.MhDraws(end-1,3);% Number of lines in the last mh files of the previous session. - LastFileNumberInThePreviousMh = sum(record.MhDraws(1:end-1,2),1);% Number of mh files in the the previous sessions. - %Test if the last mh files of the previous session were not full yet - if LastLineNumberInThePreviousMh < MAX_nruns%not full - %store starting point if whole chain needs to be redone - NewFile = ones(NumberOfBlocks,1)*LastFileNumberInThePreviousMh; - FirstLine = ones(NumberOfBlocks,1)*(LastLineNumberInThePreviousMh+1); - LastFileFullIndicator=0; - else% The last mh files of the previous session were full - %store starting point if whole chain needs to be redone - NewFile = ones(NumberOfBlocks,1)*(LastFileNumberInThePreviousMh+1); - FirstLine = ones(NumberOfBlocks,1); - LastFileFullIndicator=1; - end - else - LastLineNumberInThePreviousMh = 0; - LastFileNumberInThePreviousMh = 0; - NewFile = ones(NumberOfBlocks,1); - FirstLine = ones(NumberOfBlocks,1); - LastFileFullIndicator=1; - end - - %% Now find out what exactly needs to be redone - % 1. Check if really something needs to be done - % How many mh files should we have ? - ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1); - ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*NumberOfBlocks; - % How many mh files do we actually have ? - AllMhFiles = dir([BaseName '_mh*_blck*.mat']); - TotalNumberOfMhFiles = length(AllMhFiles); - % Quit if no crashed mcmc chain can be found as there are as many files as expected - if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles) - disp('Estimation::mcmc: It appears that you don''t need to use the mh_recover option!') - disp(' You have to edit the mod file and remove the mh_recover option') - disp(' in the estimation command') - error('Estimation::mcmc: mh_recover option not required!') - end - % 2. Something needs to be done; find out what - % Count the number of saved mh files per block. - NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1); - for b = 1:NumberOfBlocks - BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']); - NumberOfMhFilesPerBlock(b) = length(BlckMhFiles); - end - % Find FirstBlock (First block), an integer targeting the crashed mcmc chain. - FirstBlock = 1; %initialize - while FirstBlock <= NumberOfBlocks - if NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock - disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is not complete!']) - break - % The mh_recover session will start from chain FirstBlock. - else - disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is complete!']) - end - FirstBlock = FirstBlock+1; - end - - %% 3. Overwrite default settings for - % How many mh-files are saved in this block? - NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(FirstBlock); - ExistingDrawsInLastMCFile=0; %initialize: no MCMC draws of current MCMC are in file from last run - % Check whether last present file is a file included in the last MCMC run - if ~LastFileFullIndicator - if NumberOfSavedMhFilesInTheCrashedBlck==NewFile(FirstBlock) %only that last file exists, but no files from current MCMC - loaded_results=load([BaseName '_mh' int2str(NewFile(FirstBlock)) '_blck' int2str(FirstBlock) '.mat']); - %check whether that last file was filled - if size(loaded_results.x2,1)==MAX_nruns %file is full - NewFile(FirstBlock)=NewFile(FirstBlock)+1; %set first file to be created to next one - FirstLine(FirstBlock) = 1; %use first line of next file - ExistingDrawsInLastMCFile=MAX_nruns-record.MhDraws(end-1,3); - else - ExistingDrawsInLastMCFile=0; - end - end - elseif LastFileFullIndicator - ExistingDrawsInLastMCFile=0; - if NumberOfSavedMhFilesInTheCrashedBlck==NewFile(FirstBlock) %only the last file exists, but no files from current MCMC - NewFile(FirstBlock)=NewFile(FirstBlock)+1; %set first file to be created to next one - end - end - % % Correct the number of saved mh files if the crashed Metropolis was not the first session (so - % % that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain - % % of the current session). - % if OldMhExists - % NumberOfSavedMhFilesInTheCrashedBlck = NumberOfSavedMhFilesInTheCrashedBlck - LastFileNumberInThePreviousMh; - % end - % NumberOfSavedMhFiles = NumberOfSavedMhFilesInTheCrashedBlck+LastFileNumberInThePreviousMh; - - % Correct initial conditions. - if NumberOfSavedMhFilesInTheCrashedBlck