1function [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds) 2 3% function [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds) 4% initialization of posterior samplers 5% 6% INPUTS 7% posterior_sampler_options: posterior sampler options 8% options_: structure storing the options 9 10% OUTPUTS 11% posterior_sampler_options: checked posterior sampler options 12% 13% SPECIAL REQUIREMENTS 14% none 15 16% Copyright (C) 2015-2017 Dynare Team 17% 18% This file is part of Dynare. 19% 20% Dynare is free software: you can redistribute it and/or modify 21% it under the terms of the GNU General Public License as published by 22% the Free Software Foundation, either version 3 of the License, or 23% (at your option) any later version. 24% 25% Dynare is distributed in the hope that it will be useful, 26% but WITHOUT ANY WARRANTY; without even the implied warranty of 27% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 28% GNU General Public License for more details. 29% 30% You should have received a copy of the GNU General Public License 31% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 32 33 34init=0; 35if isempty(posterior_sampler_options) 36 init=1; 37end 38 39if init 40 % set default options and user defined options 41 posterior_sampler_options.posterior_sampling_method = options_.posterior_sampler_options.posterior_sampling_method; 42 posterior_sampler_options.bounds = bounds; 43 44 switch posterior_sampler_options.posterior_sampling_method 45 46 case 'random_walk_metropolis_hastings' 47 posterior_sampler_options.parallel_bar_refresh_rate=50; 48 posterior_sampler_options.serial_bar_refresh_rate=3; 49 posterior_sampler_options.parallel_bar_title='RWMH'; 50 posterior_sampler_options.serial_bar_title='RW Metropolis-Hastings'; 51 52 % default options 53 posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.rwmh); 54 55 % user defined options 56 if ~isempty(options_.posterior_sampler_options.sampling_opt) 57 options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt); 58 for i=1:rows(options_list) 59 switch options_list{i,1} 60 61 case 'proposal_distribution' 62 if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ... 63 strcmpi(options_list{i,2}, 'rand_multivariate_normal')) 64 error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ... 65 'rand_multivariate_student or rand_multivariate_normal as options']); 66 else 67 posterior_sampler_options.proposal_distribution=options_list{i,2}; 68 end 69 70 71 case 'student_degrees_of_freedom' 72 if options_list{i,2} <= 0 73 error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument'); 74 else 75 posterior_sampler_options.student_degrees_of_freedom=options_list{i,2}; 76 end 77 78 case 'use_mh_covariance_matrix' 79 % indicates to use the covariance matrix from previous iterations to 80 % define the covariance of the proposal distribution 81 % default = 0 82 posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2}; 83 options_.use_mh_covariance_matrix = options_list{i,2}; 84 case 'scale_file' 85 % load optimal_mh_scale parameter if previous run was with mode_compute=6 86 % will overwrite jscale from set_prior.m 87 if exist(options_list{i,2},'file') || exist([options_list{i,2},'.mat'],'file') 88 tmp = load(options_list{i,2},'Scale'); 89 global bayestopt_ 90 bayestopt_.mh_jscale = tmp.Scale; 91 options_.mh_jscale = tmp.Scale; 92 bayestopt_.jscale = ones(size(bounds.lb,1),1)*tmp.Scale; 93 % options_.mh_init_scale = 2*options_.mh_jscale; 94 else 95 error('initial_estimation_checks:: The specified mh_scale_file does not exist.') 96 end 97 case 'save_tmp_file' 98 posterior_sampler_options.save_tmp_file = options_list{i,2}; 99 otherwise 100 warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!']) 101 end 102 end 103 end 104 105 case 'tailored_random_block_metropolis_hastings' 106 posterior_sampler_options.parallel_bar_refresh_rate=5; 107 posterior_sampler_options.serial_bar_refresh_rate=1; 108 posterior_sampler_options.parallel_bar_title='TaRB-MH'; 109 posterior_sampler_options.serial_bar_title='TaRB Metropolis-Hastings'; 110 111 % default options 112 posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.tarb); 113 114 % user defined options 115 if ~isempty(options_.posterior_sampler_options.sampling_opt) 116 options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt); 117 for i=1:rows(options_list) 118 119 switch options_list{i,1} 120 121 case 'proposal_distribution' 122 if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ... 123 strcmpi(options_list{i,2}, 'rand_multivariate_normal')) 124 error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ... 125 'rand_multivariate_student or rand_multivariate_normal as options']); 126 else 127 posterior_sampler_options.proposal_distribution=options_list{i,2}; 128 end 129 130 131 case 'student_degrees_of_freedom' 132 if options_list{i,2} <= 0 133 error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument'); 134 else 135 posterior_sampler_options.student_degrees_of_freedom=options_list{i,2}; 136 end 137 138 case 'mode_compute' 139 posterior_sampler_options.mode_compute=options_list{i,2}; 140 141 case 'optim' 142 posterior_sampler_options.optim_opt=options_list{i,2}; 143 144 case 'new_block_probability' 145 if options_list{i,2}<0 || options_list{i,2}>1 146 error('check_posterior_sampler_options:: The tarb new_block_probability must be between 0 and 1!') 147 else 148 posterior_sampler_options.new_block_probability=options_list{i,2}; 149 end 150 case 'scale_file' 151 % load optimal_mh_scale parameter if previous run was with mode_compute=6 152 % will overwrite jscale from set_prior.m 153 if exist(options_list{i,2},'file') || exist([options_list{i,2},'.mat'],'file') 154 tmp = load(options_list{i,2},'Scale'); 155 global bayestopt_ 156 bayestopt_.mh_jscale = tmp.Scale; 157 options_.mh_jscale = tmp.Scale; 158 bayestopt_.jscale = ones(size(bounds.lb,1),1)*tmp.Scale; 159 % options_.mh_init_scale = 2*options_.mh_jscale; 160 else 161 error('initial_estimation_checks:: The specified scale_file does not exist.') 162 end 163 case 'save_tmp_file' 164 posterior_sampler_options.save_tmp_file = options_list{i,2}; 165 166 otherwise 167 warning(['tarb_sampler: Unknown option (' options_list{i,1} ')!']) 168 169 end 170 171 end 172 173 end 174 175 case 'independent_metropolis_hastings' 176 posterior_sampler_options.parallel_bar_refresh_rate=50; 177 posterior_sampler_options.serial_bar_refresh_rate=3; 178 posterior_sampler_options.parallel_bar_title='IMH'; 179 posterior_sampler_options.serial_bar_title='Ind. Metropolis-Hastings'; 180 181 % default options 182 posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.imh); 183 184 % user defined options 185 if ~isempty(options_.posterior_sampler_options.sampling_opt) 186 options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt); 187 for i=1:rows(options_list) 188 switch options_list{i,1} 189 190 case 'proposal_distribution' 191 if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ... 192 strcmpi(options_list{i,2}, 'rand_multivariate_normal')) 193 error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ... 194 'rand_multivariate_student or rand_multivariate_normal as options']); 195 else 196 posterior_sampler_options.proposal_distribution=options_list{i,2}; 197 end 198 199 200 case 'student_degrees_of_freedom' 201 if options_list{i,2} <= 0 202 error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument'); 203 else 204 posterior_sampler_options.student_degrees_of_freedom=options_list{i,2}; 205 end 206 207 case 'use_mh_covariance_matrix' 208 % indicates to use the covariance matrix from previous iterations to 209 % define the covariance of the proposal distribution 210 % default = 0 211 posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2}; 212 options_.use_mh_covariance_matrix = options_list{i,2}; 213 214 case 'save_tmp_file' 215 posterior_sampler_options.save_tmp_file = options_list{i,2}; 216 217 otherwise 218 warning(['imh_sampler: Unknown option (' options_list{i,1} ')!']) 219 end 220 end 221 end 222 223 224 case 'slice' 225 posterior_sampler_options.parallel_bar_refresh_rate=1; 226 posterior_sampler_options.serial_bar_refresh_rate=1; 227 posterior_sampler_options.parallel_bar_title='SLICE'; 228 posterior_sampler_options.serial_bar_title='SLICE'; 229 230 % default options 231 posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.slice); 232 233 % user defined options 234 if ~isempty(options_.posterior_sampler_options.sampling_opt) 235 options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt); 236 for i=1:rows(options_list) 237 switch options_list{i,1} 238 case 'rotated' 239 % triggers rotated slice iterations using a covariance 240 % matrix from initial burn-in iterations 241 % must be associated with: 242 % <use_mh_covariance_matrix> or <slice_initialize_with_mode> 243 % default = 0 244 posterior_sampler_options.rotated = options_list{i,2}; 245 246 case 'mode' 247 % for multimodal posteriors, provide the list of modes as a 248 % matrix, ordered by column, i.e. [x1 x2 x3] for three 249 % modes x1 x2 x3 250 % MR note: not sure this is possible with the 251 % read_key_value_string ??? 252 % if this is not possible <mode_files> does to job in any case 253 % This will automatically trigger <rotated> 254 % default = [] 255 tmp_mode = options_list{i,2}; 256 for j=1:size(tmp_mode,2) 257 posterior_sampler_options.mode(j).m = tmp_mode(:,j); 258 end 259 260 case 'mode_files' 261 % for multimodal posteriors provide the name of 262 % a file containing a variable array xparams = [nparam * nmodes] 263 % one column per mode. With this info, the code will automatically 264 % set the <mode> option. 265 % This will automatically trigger <rotated> 266 % default = [] 267 posterior_sampler_options.mode_files = options_list{i,2}; 268 269 case 'slice_initialize_with_mode' 270 % the default for slice is to set mode_compute = 0 in the 271 % preprocessor and start the chain(s) from a random location in the prior. 272 % This option first runs the optimizer and then starts the 273 % chain from the mode. Associated with optios <rotated>, it will 274 % use invhess from the mode to perform rotated slice 275 % iterations. 276 % default = 0 277 posterior_sampler_options.slice_initialize_with_mode = options_list{i,2}; 278 279 case 'initial_step_size' 280 % sets the initial size of the interval in the STEPPING-OUT PROCEDURE 281 % the initial_step_size must be a real number in [0, 1], 282 % and it sets the size as a proportion of the prior bounds, 283 % i.e. the size will be initial_step_size*(UB-LB) 284 % slice sampler requires prior_truncation > 0! 285 % default = 0.8 286 if options_list{i,2}<=0 || options_list{i,2}>=1 287 error('check_posterior_sampler_options:: slice initial_step_size must be between 0 and 1') 288 else 289 posterior_sampler_options.initial_step_size=options_list{i,2}; 290 end 291 case 'use_mh_covariance_matrix' 292 % in association with <rotated> indicates to use the 293 % covariance matrix from previous iterations to define the 294 % rotated slice 295 % default = 0 296 posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2}; 297 options_.use_mh_covariance_matrix = options_list{i,2}; 298 299 case 'save_tmp_file' 300 posterior_sampler_options.save_tmp_file = options_list{i,2}; 301 302 otherwise 303 warning(['slice_sampler: Unknown option (' options_list{i,1} ')!']) 304 end 305 end 306 end 307 308 % slice posterior sampler does not require mode or hessian to run 309 % needs to be set to 1 to skip parts in dynare_estimation_1.m 310 % requiring posterior maximization/calibrated smoother before MCMC 311 options_.mh_posterior_mode_estimation=1; 312 313 if ~ posterior_sampler_options.slice_initialize_with_mode 314 % by default, slice sampler should trigger 315 % mode_compute=0 and 316 % mh_replic=100 (much smaller than the default mh_replic=20000 of RWMH) 317 options_.mode_compute = 0; 318 options_.cova_compute = 0; 319 else 320 if (isequal(options_.mode_compute,0) && isempty(options_.mode_file) ) 321 skipline() 322 disp('check_posterior_sampler_options:: You have specified the option "slice_initialize_with_mode"') 323 disp('check_posterior_sampler_options:: to initialize the slice sampler using mode information') 324 disp('check_posterior_sampler_options:: but no mode file nor posterior maximization is selected,') 325 error('check_posterior_sampler_options:: The option "slice_initialize_with_mode" is inconsistent with mode_compute=0 or empty mode_file.') 326 else 327 options_.mh_posterior_mode_estimation=0; 328 end 329 end 330 331 if any(isinf(bounds.lb)) || any(isinf(bounds.ub)) 332 skipline() 333 disp('some priors are unbounded and prior_trunc is set to zero') 334 error('The option "slice" is inconsistent with prior_trunc=0.') 335 end 336 337 % moreover slice must be associated to: 338 % options_.mh_posterior_mode_estimation = 0; 339 % this is done below, but perhaps preprocessing should do this? 340 341 if ~isempty(posterior_sampler_options.mode) 342 % multimodal case 343 posterior_sampler_options.rotated = 1; 344 posterior_sampler_options.WR=[]; 345 end 346 % posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]); 347 348 349 posterior_sampler_options.W1=posterior_sampler_options.initial_step_size*(bounds.ub-bounds.lb); 350 if options_.load_mh_file 351 posterior_sampler_options.slice_initialize_with_mode = 0; 352 else 353 if ~posterior_sampler_options.slice_initialize_with_mode 354 posterior_sampler_options.invhess=[]; 355 end 356 end 357 358 if ~isempty(posterior_sampler_options.mode_files) % multimodal case 359 modes = posterior_sampler_options.mode_files; % these can be also mean files from previous parallel slice chains 360 load(modes, 'xparams') 361 if size(xparams,2)<2 362 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.']) 363 end 364 for j=1:size(xparams,2) 365 mode(j).m=xparams(:,j); 366 end 367 posterior_sampler_options.mode = mode; 368 posterior_sampler_options.rotated = 1; 369 posterior_sampler_options.WR=[]; 370 end 371 372 otherwise 373 error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method); 374 end 375 376 return 377end 378 379% here are all samplers requiring a proposal distribution 380if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice') 381 if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix) 382 if strcmp('hessian',options_.MCMC_jumping_covariance) 383 skipline() 384 disp('check_posterior_sampler_options:: I cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed') 385 disp('check_posterior_sampler_options:: or there is no previous MCMC to load ') 386 error('check_posterior_sampler_options:: MCMC cannot start') 387 end 388 end 389end 390 391if options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix 392 [~, invhess] = compute_mh_covariance_matrix; 393 posterior_sampler_options.invhess = invhess; 394end 395 396 397 398% check specific options for slice sampler 399if strcmp(posterior_sampler_options.posterior_sampling_method,'slice') 400 invhess = posterior_sampler_options.invhess; 401 if posterior_sampler_options.rotated 402 if isempty(posterior_sampler_options.mode_files) && isempty(posterior_sampler_options.mode) % rotated unimodal 403 if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix) 404 skipline() 405 disp('check_posterior_sampler_options:: I cannot start rotated slice sampler because') 406 disp('check_posterior_sampler_options:: there is no previous MCMC to load ') 407 disp('check_posterior_sampler_options:: or the Hessian at the mode is not computed.') 408 error('check_posterior_sampler_options:: Rotated slice cannot start') 409 end 410 if isempty(invhess) 411 error('check_posterior_sampler_options:: This error should not occur, please contact developers.') 412 end 413 % % % if options_.load_mh_file && options_.use_mh_covariance_matrix, 414 % % % [~, invhess] = compute_mh_covariance_matrix; 415 % % % posterior_sampler_options.invhess = invhess; 416 % % % end 417 [V1, D]=eig(invhess); 418 posterior_sampler_options.V1=V1; 419 posterior_sampler_options.WR=sqrt(diag(D))*3; 420 end 421 else 422 if ~options_.load_mh_file && ~posterior_sampler_options.slice_initialize_with_mode 423 posterior_sampler_options.invhess=[]; 424 end 425 end 426 % needs to be re-set to zero otherwise posterior analysis is filtered 427 % out in dynare_estimation_1.m 428 options_.mh_posterior_mode_estimation = 0; 429end 430 431return 432 433function posterior_sampler_options = add_fields_(posterior_sampler_options, sampler_options) 434 435fnam = fieldnames(sampler_options); 436for j=1:length(fnam) 437 posterior_sampler_options.(fnam{j}) = sampler_options.(fnam{j}); 438end 439