Select Git revision
get_error_message.m
check_posterior_sampler_options.m 21.03 KiB
function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds, bayestopt_)
% function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds, bayestopt_)
% initialization of posterior samplers
%
% INPUTS
% posterior_sampler_options: posterior sampler options
% options_: structure storing the options
% bounds: structure containing prior bounds
% bayestopt_: structure storing information about priors
% OUTPUTS
% posterior_sampler_options: checked posterior sampler options
% options_: structure storing the options
% bayestopt_: structure storing information about priors
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015-2022 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 <https://www.gnu.org/licenses/>.
init=0;
if isempty(posterior_sampler_options)
init=1;
end
if init
% set default options and user defined options
posterior_sampler_options.posterior_sampling_method = options_.posterior_sampler_options.posterior_sampling_method;
posterior_sampler_options.bounds = bounds;
switch posterior_sampler_options.posterior_sampling_method
case 'random_walk_metropolis_hastings'
posterior_sampler_options.parallel_bar_refresh_rate=50;
posterior_sampler_options.serial_bar_refresh_rate=3;
posterior_sampler_options.parallel_bar_title='RWMH';
posterior_sampler_options.serial_bar_title='RW Metropolis-Hastings';
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.rwmh);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'proposal_distribution'
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
'rand_multivariate_student or rand_multivariate_normal as options']);
else posterior_sampler_options.proposal_distribution=options_list{i,2};
end
case 'student_degrees_of_freedom'
if options_list{i,2} <= 0
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
else
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end
case 'use_mh_covariance_matrix'
% indicates to use the covariance matrix from previous iterations to
% define the covariance of the proposal distribution
% default = 0
posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2};
options_.use_mh_covariance_matrix = options_list{i,2};
case 'scale_file'
% load optimal_mh_scale parameter if previous run was with mode_compute=6
% will overwrite jscale from set_prior.m
if exist(options_list{i,2},'file') || exist([options_list{i,2},'.mat'],'file')
tmp = load(options_list{i,2},'Scale');
bayestopt_.mh_jscale = tmp.Scale;
options_.mh_jscale = tmp.Scale;
bayestopt_.jscale = ones(size(bounds.lb,1),1)*tmp.Scale;
% options_.mh_init_scale = 2*options_.mh_jscale;
else
error('initial_estimation_checks:: The specified mh_scale_file does not exist.')
end
case 'save_tmp_file'
posterior_sampler_options.save_tmp_file = options_list{i,2};
otherwise
warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
end
end
end
case 'tailored_random_block_metropolis_hastings'
posterior_sampler_options.parallel_bar_refresh_rate=5;
posterior_sampler_options.serial_bar_refresh_rate=1;
posterior_sampler_options.parallel_bar_title='TaRB-MH';
posterior_sampler_options.serial_bar_title='TaRB Metropolis-Hastings';
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.tarb);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'proposal_distribution'
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
'rand_multivariate_student or rand_multivariate_normal as options']);
else
posterior_sampler_options.proposal_distribution=options_list{i,2};
end
case 'student_degrees_of_freedom'
if options_list{i,2} <= 0
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
else
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end
case 'mode_compute'
posterior_sampler_options.mode_compute=options_list{i,2};
case 'optim'
posterior_sampler_options.optim_opt=options_list{i,2};
case 'new_block_probability'
if options_list{i,2}<0 || options_list{i,2}>1
error('check_posterior_sampler_options:: The tarb new_block_probability must be between 0 and 1!')
else
posterior_sampler_options.new_block_probability=options_list{i,2};
end
case 'scale_file'
% load optimal_mh_scale parameter if previous run was with mode_compute=6
% will overwrite jscale from set_prior.m
if exist(options_list{i,2},'file') || exist([options_list{i,2},'.mat'],'file')
tmp = load(options_list{i,2},'Scale');
bayestopt_.mh_jscale = tmp.Scale;
options_.mh_jscale = tmp.Scale;
bayestopt_.jscale = ones(size(bounds.lb,1),1)*tmp.Scale;
% options_.mh_init_scale = 2*options_.mh_jscale;
else
error('initial_estimation_checks:: The specified scale_file does not exist.')
end
case 'save_tmp_file'
posterior_sampler_options.save_tmp_file = options_list{i,2};
otherwise
warning(['tarb_sampler: Unknown option (' options_list{i,1} ')!'])
end
end
end
case 'independent_metropolis_hastings'
posterior_sampler_options.parallel_bar_refresh_rate=50;
posterior_sampler_options.serial_bar_refresh_rate=3;
posterior_sampler_options.parallel_bar_title='IMH';
posterior_sampler_options.serial_bar_title='Ind. Metropolis-Hastings';
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.imh);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'proposal_distribution'
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
'rand_multivariate_student or rand_multivariate_normal as options']);
else
posterior_sampler_options.proposal_distribution=options_list{i,2};
end
case 'student_degrees_of_freedom'
if options_list{i,2} <= 0
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
else
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end
case 'use_mh_covariance_matrix'
% indicates to use the covariance matrix from previous iterations to % define the covariance of the proposal distribution
% default = 0
posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2};
options_.use_mh_covariance_matrix = options_list{i,2};
case 'save_tmp_file'
posterior_sampler_options.save_tmp_file = options_list{i,2};
otherwise
warning(['imh_sampler: Unknown option (' options_list{i,1} ')!'])
end
end
end
case 'slice'
posterior_sampler_options.parallel_bar_refresh_rate=1;
posterior_sampler_options.serial_bar_refresh_rate=1;
posterior_sampler_options.parallel_bar_title='SLICE';
posterior_sampler_options.serial_bar_title='SLICE';
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.slice);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'rotated'
% triggers rotated slice iterations using a covariance
% matrix from initial burn-in iterations
% must be associated with:
% <use_mh_covariance_matrix> or <slice_initialize_with_mode>
% default = 0
posterior_sampler_options.rotated = options_list{i,2};
case 'mode'
% for multimodal posteriors, provide the list of modes as a
% matrix, ordered by column, i.e. [x1 x2 x3] for three
% modes x1 x2 x3
% MR note: not sure this is possible with the
% read_key_value_string ???
% if this is not possible <mode_files> does to job in any case
% This will automatically trigger <rotated>
% default = []
tmp_mode = options_list{i,2};
for j=1:size(tmp_mode,2)
posterior_sampler_options.mode(j).m = tmp_mode(:,j);
end
case 'mode_files'
% for multimodal posteriors provide the name of
% a file containing a variable array xparams = [nparam * nmodes]
% one column per mode. With this info, the code will automatically
% set the <mode> option.
% This will automatically trigger <rotated>
% default = []
posterior_sampler_options.mode_files = options_list{i,2};
case 'slice_initialize_with_mode'
% the default for slice is to set mode_compute = 0 in the
% preprocessor and start the chain(s) from a random location in the prior.
% This option first runs the optimizer and then starts the
% chain from the mode. Associated with optios <rotated>, it will
% use invhess from the mode to perform rotated slice
% iterations.
% default = 0
posterior_sampler_options.slice_initialize_with_mode = options_list{i,2};
case 'initial_step_size'
% sets the initial size of the interval in the STEPPING-OUT PROCEDURE
% the initial_step_size must be a real number in [0, 1],
% and it sets the size as a proportion of the prior bounds,
% i.e. the size will be initial_step_size*(UB-LB)
% slice sampler requires prior_truncation > 0!
% default = 0.8
if options_list{i,2}<=0 || options_list{i,2}>=1
error('check_posterior_sampler_options:: slice initial_step_size must be between 0 and 1')
else
posterior_sampler_options.initial_step_size=options_list{i,2};
end
case 'use_mh_covariance_matrix'
% in association with <rotated> indicates to use the
% covariance matrix from previous iterations to define the
% rotated slice
% default = 0
posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2};
options_.use_mh_covariance_matrix = options_list{i,2};
case 'save_tmp_file'
posterior_sampler_options.save_tmp_file = options_list{i,2};
otherwise
warning(['slice_sampler: Unknown option (' options_list{i,1} ')!'])
end
end
end
% slice posterior sampler does not require mode or hessian to run
% needs to be set to 1 to skip parts in dynare_estimation_1.m
% requiring posterior maximization/calibrated smoother before MCMC
options_.mh_posterior_mode_estimation=true;
if ~ posterior_sampler_options.slice_initialize_with_mode
% by default, slice sampler should trigger
% mode_compute=0 and
% mh_replic=100 (much smaller than the default mh_replic=20000 of RWMH)
options_.mode_compute = 0;
options_.cova_compute = 0;
else
if (isequal(options_.mode_compute,0) && isempty(options_.mode_file) )
skipline()
disp('check_posterior_sampler_options:: You have specified the option "slice_initialize_with_mode"')
disp('check_posterior_sampler_options:: to initialize the slice sampler using mode information')
disp('check_posterior_sampler_options:: but no mode file nor posterior maximization is selected,')
error('check_posterior_sampler_options:: The option "slice_initialize_with_mode" is inconsistent with mode_compute=0 or empty mode_file.')
else
options_.mh_posterior_mode_estimation=false;
end
end
if any(isinf(bounds.lb)) || any(isinf(bounds.ub))
skipline()
disp('some priors are unbounded and prior_trunc is set to zero')
error('The option "slice" is inconsistent with prior_trunc=0.')
end
% moreover slice must be associated to:
% options_.mh_posterior_mode_estimation = false;
% this is done below, but perhaps preprocessing should do this?
if ~isempty(posterior_sampler_options.mode)
% multimodal case
posterior_sampler_options.rotated = 1;
posterior_sampler_options.WR=[];
end
% posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]);
posterior_sampler_options.W1=posterior_sampler_options.initial_step_size*(bounds.ub-bounds.lb);
if options_.load_mh_file
posterior_sampler_options.slice_initialize_with_mode = 0;
else
if ~posterior_sampler_options.slice_initialize_with_mode
posterior_sampler_options.invhess=[];
end
end
if ~isempty(posterior_sampler_options.mode_files) % multimodal case
modes = posterior_sampler_options.mode_files; % these can be also mean files from previous parallel slice chains
load(modes, 'xparams')
if size(xparams,2)<2
error(['check_posterior_sampler_options:: Variable xparams loaded in file <' modes '> has size [' int2str(size(xparams,1)) 'x' int2str(size(xparams,2)) ']: it must contain at least two columns, to allow multi-modal sampling.'])
end
for j=1:size(xparams,2)
mode(j).m=xparams(:,j);
end
posterior_sampler_options.mode = mode;
posterior_sampler_options.rotated = 1;
posterior_sampler_options.WR=[];
end
otherwise
error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
end
return
end
% here are all samplers requiring a proposal distribution
if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix)
if strcmp('hessian',options_.MCMC_jumping_covariance)
skipline()
disp('check_posterior_sampler_options:: I cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed')
disp('check_posterior_sampler_options:: or there is no previous MCMC to load ')
error('check_posterior_sampler_options:: MCMC cannot start')
end
end
end
if options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix
[~, invhess] = compute_mh_covariance_matrix;
posterior_sampler_options.invhess = invhess;
end
% check specific options for slice sampler
if strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
invhess = posterior_sampler_options.invhess;
if posterior_sampler_options.rotated
if isempty(posterior_sampler_options.mode_files) && isempty(posterior_sampler_options.mode) % rotated unimodal
if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix)
skipline()
disp('check_posterior_sampler_options:: I cannot start rotated slice sampler because')
disp('check_posterior_sampler_options:: there is no previous MCMC to load ')
disp('check_posterior_sampler_options:: or the Hessian at the mode is not computed.')
error('check_posterior_sampler_options:: Rotated slice cannot start')
end
if isempty(invhess)
error('check_posterior_sampler_options:: This error should not occur, please contact developers.')
end
% % % if options_.load_mh_file && options_.use_mh_covariance_matrix,
% % % [~, invhess] = compute_mh_covariance_matrix;
% % % posterior_sampler_options.invhess = invhess;
% % % end
[V1, D]=eig(invhess);
posterior_sampler_options.V1=V1; posterior_sampler_options.WR=sqrt(diag(D))*3;
end
else
if ~options_.load_mh_file && ~posterior_sampler_options.slice_initialize_with_mode
posterior_sampler_options.invhess=[];
end
end
% needs to be re-set to zero otherwise posterior analysis is filtered
% out in dynare_estimation_1.m
options_.mh_posterior_mode_estimation = false;
end
return
function posterior_sampler_options = add_fields_(posterior_sampler_options, sampler_options)
fnam = fieldnames(sampler_options);
for j=1:length(fnam)
posterior_sampler_options.(fnam{j}) = sampler_options.(fnam{j});
end