2016-05-13 21:35:59 +02:00
function [posterior_sampler_options, options_] = check_posterior_sampler_options ( posterior_sampler_options, options_, bounds)
% function [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds)
% initialization of posterior samplers
%
% INPUTS
% posterior_sampler_options: posterior sampler options
% options_: structure storing the options
% OUTPUTS
% posterior_sampler_options: checked posterior sampler options
%
% SPECIAL REQUIREMENTS
% none
2017-05-16 14:11:15 +02:00
% Copyright (C) 2015-2017 Dynare Team
2016-05-13 21:35:59 +02:00
%
% 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 <http://www.gnu.org/licenses/>.
init = 0 ;
2017-05-16 12:42:01 +02:00
if isempty ( posterior_sampler_options )
2016-05-13 21:35:59 +02:00
init = 1 ;
end
2017-05-16 12:42:01 +02:00
if init
2016-05-13 21:35:59 +02:00
% 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 ;
2017-05-16 15:10:20 +02:00
2016-05-13 21:35:59 +02:00
switch posterior_sampler_options . posterior_sampling_method
2017-05-16 15:10:20 +02:00
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' ) ;
global bayestopt_
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.' )
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
case ' save_tmp_file'
posterior_sampler_options . save_tmp_file = options_list { i , 2 } ;
otherwise
warning ( [ ' rwmh_sampler: Unknown option (' options_list { i , 1 } ' )!' ] )
2016-05-13 21:35:59 +02:00
end
end
2017-05-16 15:10:20 +02:00
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 } ;
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
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' ) ;
global bayestopt_
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 } ' )!' ] )
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
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 } ;
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
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 } ' )!' ] )
2016-05-13 21:35:59 +02:00
end
end
2017-05-16 15:10:20 +02:00
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 } ;
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
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 } ' )!' ] )
2016-05-13 21:35:59 +02:00
end
end
2017-05-16 15:10:20 +02:00
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 = 1 ;
if ~ posterior_sampler_options . slice_initialize_with_mode
2016-05-13 21:35:59 +02:00
% by default, slice sampler should trigger
% mode_compute=0 and
% mh_replic=100 (much smaller than the default mh_replic=20000 of RWMH)
2017-05-16 15:10:20 +02:00
options_ . mode_compute = 0 ;
options_ . cova_compute = 0 ;
else
if ( isequal ( options_ . mode_compute , 0 ) && isempty ( options_ . mode_file ) )
2016-05-13 21:35:59 +02:00
skipline ( )
2017-05-16 15:10:20 +02:00
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 = 0 ;
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
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 = 0;
% 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 = [ ] ;
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
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.' ] )
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
for j = 1 : size ( xparams , 2 )
mode ( j ) . m = xparams ( : , j ) ;
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
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 ) ;
2016-05-13 21:35:59 +02:00
end
2017-05-16 15:10:20 +02:00
2016-05-13 21:35:59 +02:00
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 )
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
2017-05-16 12:42:01 +02:00
if options_ . load_mh_file && posterior_sampler_options . use_mh_covariance_matrix
2018-11-13 17:58:42 +01:00
[ ~ , invhess ] = compute_mh_covariance_matrix ;
2016-05-13 21:35:59 +02:00
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 ;
2017-05-16 12:42:01 +02:00
if posterior_sampler_options . rotated
if isempty ( posterior_sampler_options . mode_files ) && isempty ( posterior_sampler_options . mode ) % rotated unimodal
2016-05-13 21:35:59 +02:00
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,
2018-11-13 17:58:42 +01:00
% % % [~, invhess] = compute_mh_covariance_matrix;
2016-05-13 21:35:59 +02:00
% % % posterior_sampler_options.invhess = invhess;
% % % end
2016-10-04 14:56:41 +02:00
[ V1 , D ] = eig ( invhess ) ;
2016-05-13 21:35:59 +02:00
posterior_sampler_options . V1 = V1 ;
posterior_sampler_options . WR = sqrt ( diag ( D ) ) * 3 ;
end
2017-03-22 23:06:19 +01:00
else
2017-05-16 12:42:01 +02:00
if ~ options_ . load_mh_file && ~ posterior_sampler_options . slice_initialize_with_mode
2017-03-22 23:06:19 +01:00
posterior_sampler_options . invhess = [ ] ;
end
2016-05-13 21:35:59 +02:00
end
% needs to be re-set to zero otherwise posterior analysis is filtered
% out in dynare_estimation_1.m
options_ . mh_posterior_mode_estimation = 0 ;
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