From c30b6f20a62b902974de6546931a7e89da745da7 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 25 Sep 2015 12:54:52 +0200 Subject: [PATCH] Added routines that factor out posterior samplers into the iteration routine posterior_sampler_iteration.m. It generalizes the random_walk_metropolis_hastings.m structure but adds two new items: - check_posterior_sampler_options.m which is a wrapper for method specific options (stored in options_.posterior_options) - posterior_sampler_iteration.m which is a wrapper of the individual iteration of any posterior sampler. It already embeds, with parallelization implementation: - random_walk metropolis hastings - independent metropolis hastings - slice sampler It also contains a quick fix for TARB, although the latter should be embedded in posterior_sampler_iteration.m as well. Any new posterior sampler can be simply added to posterior_sampler_iteration.m, with wrapper to specific iterations in check_posterior_sampler_options.m --- matlab/check_posterior_sampler_options.m | 116 +++++++ matlab/posterior_sampler.m | 196 +++++++++++ matlab/posterior_sampler_core.m | 241 ++++++++++++++ matlab/posterior_sampler_initialization.m | 384 ++++++++++++++++++++++ matlab/posterior_sampler_iteration.m | 93 ++++++ 5 files changed, 1030 insertions(+) create mode 100644 matlab/check_posterior_sampler_options.m create mode 100644 matlab/posterior_sampler.m create mode 100644 matlab/posterior_sampler_core.m create mode 100644 matlab/posterior_sampler_initialization.m create mode 100644 matlab/posterior_sampler_iteration.m diff --git a/matlab/check_posterior_sampler_options.m b/matlab/check_posterior_sampler_options.m new file mode 100644 index 000000000..5f0cd7b7d --- /dev/null +++ b/matlab/check_posterior_sampler_options.m @@ -0,0 +1,116 @@ +function [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_) + +% function posterior_sampler_options = check_posterior_sampler_options(posterior_sampler_options, options_) +% 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 + +% Copyright (C) 2015 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 . + +posterior_sampler_options.posterior_sampling_method = options_.posterior_sampling_method; +posterior_sampler_options.proposal_distribution = options_.proposal_distribution; +bounds = posterior_sampler_options.bounds; +invhess = posterior_sampler_options.invhess; + +% here are all samplers requiring a proposal distribution +if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice') + if ~options_.cova_compute + error('I Cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed.') + end + if options_.load_mh_file && options_.use_mh_covariance_matrix, + invhess = compute_mh_covariance_matrix; + posterior_sampler_options.invhess = invhess; + end + posterior_sampler_options.parallel_bar_refresh_rate=50; + posterior_sampler_options.serial_bar_refresh_rate=3; + posterior_sampler_options.parallel_bar_title='MH'; + posterior_sampler_options.serial_bar_title='Metropolis-Hastings'; + posterior_sampler_options.save_tmp_file=1; +end + + +% check specific options for slice sampler +if strcmp(options_.posterior_sampling_method,'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'; + posterior_sampler_options.save_tmp_file=1; + posterior_sampler_options = set_default_option(posterior_sampler_options,'rotated',0); + posterior_sampler_options = set_default_option(posterior_sampler_options,'slice_initialize_with_mode',0); + posterior_sampler_options = set_default_option(posterior_sampler_options,'use_slice_covariance_matrix',0); + posterior_sampler_options = set_default_option(posterior_sampler_options,'WR',[]); + if ~isfield(posterior_sampler_options,'mode'), + posterior_sampler_options.mode = []; + else % multimodal case + posterior_sampler_options.rotated = 1; + end + posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]); + posterior_sampler_options = set_default_option(posterior_sampler_options,'W1',0.8*(bounds.ub-bounds.lb)); + + if options_.load_mh_file, + posterior_sampler_options.slice_initialize_with_mode = 0; + end + + if posterior_sampler_options.rotated, + if isempty(posterior_sampler_options.mode_files) && ~isfield(posterior_sampler_options,'mode'), % rotated unimodal + + if ~options_.cova_compute && ~(options_.load_mh_file && options_.use_mh_covariance_matrix) + skipline() + disp('I cannot start rotated slice sampler because') + disp('there is no previous MCMC to load ') + disp('or the Hessian at the mode is not computed.') + error('Rotated slice cannot start') + end + if isempty(invhess) + error('oops! This error should not occur, please contact developers.') + end + if options_.load_mh_file && posterior_sampler_options.use_slice_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 + 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 + for j=1:length( modes ), + load(modes{j}, 'xparam1') + mode(j).m=xparam1; + end + posterior_sampler_options.mode = mode; + posterior_sampler_options.rotated = 1; + posterior_sampler_options.WR=[]; + end + options_.mh_posterior_mode_estimation = 0; + +end + + + diff --git a/matlab/posterior_sampler.m b/matlab/posterior_sampler.m new file mode 100644 index 000000000..4ec8aa140 --- /dev/null +++ b/matlab/posterior_sampler.m @@ -0,0 +1,196 @@ +function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) +% function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) +% Random Walk Metropolis-Hastings algorithm. +% +% INPUTS +% o TargetFun [char] string specifying the name of the objective +% function (posterior kernel). +% o ProposalFun [char] string specifying the name of the proposal +% density +% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values). +% o sampler_options structure +% .invhess [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 +% +% ALGORITHM +% Random-Walk Metropolis-Hastings. +% +% SPECIAL REQUIREMENTS +% None. +% +% PARALLEL CONTEXT +% The most computationally intensive part of this function may be executed +% in parallel. The code suitable to be executed in +% parallel on multi core or cluster machine (in general a 'for' cycle) +% has been removed from this function and been placed in the random_walk_metropolis_hastings_core.m funtion. +% +% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in +% parallel and called name_function.m and name_function_core.m and ii) a second set of functions used +% to manage the parallel computations. +% +% This function was the first function to be parallelized. Later, other +% functions have been parallelized using the same methodology. +% Then the comments write here can be used for all the other pairs of +% parallel functions and also for management functions. + +% Copyright (C) 2006-2015 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 . + + +% In Metropolis, we set penalty to Inf so as to reject all parameter sets triggering an error during target density computation +global objective_function_penalty_base +objective_function_penalty_base = Inf; + +vv = sampler_options.invhess; +% Initialization of the random walk metropolis-hastings chains. +[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... + posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); + +InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2); + +% Load last mh history file +load_last_mh_history_file(MetropolisFolder, ModelName); + +% Only for test parallel results!!! + +% To check the equivalence between parallel and serial computation! +% First run in serial mode, and then comment the follow line. +% save('recordSerial.mat','-struct', 'record'); + +% For parallel runs after serial runs with the abobe line active. +% TempRecord=load('recordSerial.mat'); +% record.Seeds=TempRecord.Seeds; + + + +% Snapshot of the current state of computing. It necessary for the parallel +% execution (i.e. to execute in a corretct way a portion of code remotely or +% on many cores). The mandatory variables for local/remote parallel +% computing are stored in the localVars struct. + +if options_.TaRB.use_TaRB + options_.silent_optimizer=1; %locally set optimizer to silent mode +end + +localVars = struct('TargetFun', TargetFun, ... + 'ProposalFun', ProposalFun, ... + 'xparam1', xparam1, ... + 'vv', vv, ... + 'sampler_options', sampler_options, ... + 'mh_bounds', mh_bounds, ... + 'ix2', ix2, ... + 'ilogpo2', ilogpo2, ... + 'ModelName', ModelName, ... + 'fline', fline, ... + 'npar', npar, ... + 'nruns', nruns, ... + 'NewFile', NewFile, ... + 'MAX_nruns', MAX_nruns, ... + 'd', d, ... + 'InitSizeArray',InitSizeArray, ... + 'record', record, ... + 'dataset_', dataset_, ... + 'dataset_info', dataset_info, ... + 'options_', options_, ... + 'M_',M_, ... + 'bayestopt_', bayestopt_, ... + 'estim_params_', estim_params_, ... + 'oo_', oo_,... + 'varargin',[]); + + +% User doesn't want to use parallel computing, or wants to compute a +% single chain compute Random walk Metropolis-Hastings algorithm sequentially. + +if isnumeric(options_.parallel) || (nblck-fblck)==0, + if options_.TaRB.use_TaRB + fout = TaRB_metropolis_hastings_core(localVars, fblck, nblck, 0); + else + fout = posterior_sampler_core(localVars, fblck, nblck, 0); + end + record = fout.record; + % Parallel in Local or remote machine. +else + % Global variables for parallel routines. + globalVars = struct(); + % which files have to be copied to run remotely + NamFileInput(1,:) = {'',[ModelName '_static.m']}; + NamFileInput(2,:) = {'',[ModelName '_dynamic.m']}; + if options_.steadystate_flag, + NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_steadystate.m']}; + end + if (options_.load_mh_file~=0) && any(fline>1) , + NamFileInput(length(NamFileInput)+1,:)={[M_.dname '/metropolis/'],[ModelName '_mh' int2str(NewFile(1)) '_blck*.mat']}; + end + if exist([ModelName '_optimal_mh_scale_parameter.mat']) + NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_optimal_mh_scale_parameter.mat']}; + end + % from where to get back results + % NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'}; + if options_.TaRB.use_TaRB + [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'TaRB_metropolis_hastings_core', localVars, globalVars, options_.parallel_info); + else + [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'posterior_sampler_core', localVars, globalVars, options_.parallel_info); + end + for j=1:totCPU, + offset = sum(nBlockPerCPU(1:j-1))+fblck-1; + record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j))); + record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:)=fout(j).record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:); + record.AcceptanceRatio(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptanceRatio(offset+1:sum(nBlockPerCPU(1:j))); + record.FunctionEvalPerIteration(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.FunctionEvalPerIteration(offset+1:sum(nBlockPerCPU(1:j))); + record.LastSeeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastSeeds(offset+1:sum(nBlockPerCPU(1:j))); + end + +end + +irun = fout(1).irun; +NewFile = fout(1).NewFile; + +record.MCMCConcludedSuccessfully = 1; %set indicator for successful run + +update_last_mh_history_file(MetropolisFolder, ModelName, record); + +% Provide diagnostic output +skipline() +disp(['Estimation::mcmc: Number of mh files: ' int2str(NewFile(1)) ' per block.']) +disp(['Estimation::mcmc: Total number of generated files: ' int2str(NewFile(1)*nblck) '.']) +disp(['Estimation::mcmc: Total number of iterations: ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.']) +disp(['Estimation::mcmc: Current acceptance ratio per chain: ']) +for i=1:nblck + if i<10 + disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%']) + else + disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%']) + end +end +if max(record.FunctionEvalPerIteration)>1 + disp(['Estimation::mcmc: Current function evaluations per iteration: ']) + for i=1:nblck + if i<10 + disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))]) + else + disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))]) + end + end +end \ No newline at end of file diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m new file mode 100644 index 000000000..d8a0dec71 --- /dev/null +++ b/matlab/posterior_sampler_core.m @@ -0,0 +1,241 @@ +function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatlab) +% function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatlab) +% Contains the most computationally intensive portion of code in +% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in that 'for' +% cycle are completely independent to be suitable for parallel execution. +% +% INPUTS +% o myimput [struc] The mandatory variables for local/remote +% parallel computing obtained from random_walk_metropolis_hastings.m +% function. +% o fblck and nblck [integer] The Metropolis-Hastings chains. +% o whoiam [integer] In concurrent programming a modality to refer to the different threads running in parallel is needed. +% The integer whoaim is the integer that +% allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the +% cluster. +% o ThisMatlab [integer] Allows us to distinguish between the +% 'main' Matlab, the slave Matlab worker, local Matlab, remote Matlab, +% ... Then it is the index number of this slave machine in the cluster. +% OUTPUTS +% o myoutput [struc] +% If executed without parallel, this is the original output of 'for b = +% fblck:nblck'. Otherwise, it's a portion of it computed on a specific core or +% remote machine. In this case: +% record; +% irun; +% NewFile; +% OutputFileName +% +% ALGORITHM +% Portion of Metropolis-Hastings. +% +% SPECIAL REQUIREMENTS. +% None. +% +% PARALLEL CONTEXT +% See the comments in the random_walk_metropolis_hastings.m funtion. + + +% Copyright (C) 2006-2015 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 . + +if nargin<4, + whoiam=0; +end + +% reshape 'myinputs' for local computation. +% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by: + +TargetFun=myinputs.TargetFun; +ProposalFun=myinputs.ProposalFun; +xparam1=myinputs.xparam1; +mh_bounds=myinputs.mh_bounds; +last_draw=myinputs.ix2; +last_posterior=myinputs.ilogpo2; +fline=myinputs.fline; +npar=myinputs.npar; +nruns=myinputs.nruns; +NewFile=myinputs.NewFile; +MAX_nruns=myinputs.MAX_nruns; +sampler_options=myinputs.sampler_options; +d=myinputs.d; +InitSizeArray=myinputs.InitSizeArray; +record=myinputs.record; +dataset_ = myinputs.dataset_; +dataset_info = myinputs.dataset_info; +bayestopt_ = myinputs.bayestopt_; +estim_params_ = myinputs.estim_params_; +options_ = myinputs.options_; +M_ = myinputs.M_; +oo_ = myinputs.oo_; +% Necessary only for remote computing! +if whoiam + % initialize persistent variables in priordens() + priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1); +end + +MetropolisFolder = CheckPath('metropolis',M_.dname); +ModelName = M_.fname; +BaseName = [MetropolisFolder filesep ModelName]; +save_tmp_file = sampler_options.save_tmp_file; + +options_.lik_algo = 1; +OpenOldFile = ones(nblck,1); +if strcmpi(ProposalFun,'rand_multivariate_normal') + sampler_options.n = npar; + sampler_options.ProposalDensity = 'multivariate_normal_pdf'; +elseif strcmpi(ProposalFun,'rand_multivariate_student') + sampler_options.n = options_.student_degrees_of_freedom; + sampler_options.ProposalDensity = 'multivariate_student_pdf'; +end + +% +% Now I run the (nblck-fblck+1) MCMC chains +% + +sampler_options.xparam1 = xparam1; +sampler_options.proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale); + +block_iter=0; + +for curr_block = fblck:nblck, + block_iter=block_iter+1; + try + % This will not work if the master uses a random number generator not + % available in the slave (different Matlab version or + % Matlab/Octave cluster). Therefore the trap. + % + % Set the random number generator type (the seed is useless but needed by the function) + set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed); + % Set the state of the RNG + set_dynare_random_generator_state(record.InitialSeeds(curr_block).Unifor, record.InitialSeeds(curr_block).Normal); + catch + % If the state set by master is incompatible with the slave, we only reseed + set_dynare_seed(options_.DynareRandomStreams.seed+curr_block); + end + 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); + logpo2 = zeros(InitSizeArray(curr_block),1); + end + %Prepare waiting bars + if whoiam + refresh_rate = sampler_options.parallel_bar_refresh_rate; + bar_title = sampler_options.parallel_bar_title; + prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode); + hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']); + else + refresh_rate = sampler_options.serial_bar_refresh_rate; + bar_title = sampler_options.serial_bar_title; + hh = dyn_waitbar(0,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']); + set(hh,'Name',bar_title); + end + accepted_draws_this_chain = 0; + accepted_draws_this_file = 0; + feval_this_chain = 0; + feval_this_file = 0; + draw_index_current_file = fline(curr_block); %get location of first draw in current block + draw_iter = 1; + + while draw_iter <= nruns(curr_block) + + [par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun, last_draw(curr_block,:), last_posterior(curr_block), sampler_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_); + + x2(draw_index_current_file,:) = par; + last_draw(curr_block,:) = par; + logpo2(draw_index_current_file) = logpost; + last_posterior(curr_block) = logpost; + feval_this_chain = feval_this_chain + sum(neval); + feval_this_file = feval_this_file + sum(neval); + accepted_draws_this_chain = accepted_draws_this_chain + accepted; + accepted_draws_this_file = accepted_draws_this_file + accepted; + + prtfrc = draw_iter/nruns(curr_block); + if mod(draw_iter, refresh_rate)==0 + if accepted_draws_this_chain/draw_iter==1 && sum(neval)>1 + dyn_waitbar(prtfrc,hh,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Function eval per draw %4.3f', feval_this_chain/draw_iter)]); + else + dyn_waitbar(prtfrc,hh,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/draw_iter)]); + end + if save_tmp_file + save([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'x2','logpo2'); + end + end + if (draw_index_current_file == InitSizeArray(curr_block)) || (draw_iter == nruns(curr_block)) % Now I save the simulations, either because the current file is full or the chain is done + if save_tmp_file, + delete([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']); + end + save([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'],'x2','logpo2'); + fidlog = fopen([MetropolisFolder '/metropolis.log'],'a'); + fprintf(fidlog,['\n']); + fprintf(fidlog,['%% Mh' int2str(NewFile(curr_block)) 'Blck' int2str(curr_block) ' (' datestr(now,0) ')\n']); + fprintf(fidlog,' \n'); + fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']); + fprintf(fidlog,[' Acceptance ratio......: ' num2str(accepted_draws_this_file/length(logpo2)) '\n']); + fprintf(fidlog,[' Feval per iteration...: ' num2str(feval_this_file/length(logpo2)) '\n']); + fprintf(fidlog,[' Posterior mean........:\n']); + for i=1:length(x2(1,:)) + fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']); + end + fprintf(fidlog,[' log2po:' num2str(mean(logpo2)) '\n']); + fprintf(fidlog,[' Minimum value.........:\n']); + for i=1:length(x2(1,:)) + fprintf(fidlog,[' params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']); + end + fprintf(fidlog,[' log2po:' num2str(min(logpo2)) '\n']); + fprintf(fidlog,[' Maximum value.........:\n']); + for i=1:length(x2(1,:)) + fprintf(fidlog,[' params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']); + end + fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']); + fprintf(fidlog,' \n'); + fclose(fidlog); + accepted_draws_this_file = 0; + feval_this_file = 0; + if draw_iter == nruns(curr_block) % I record the last draw... + record.LastParameters(curr_block,:) = x2(end,:); + record.LastLogPost(curr_block) = logpo2(end); + end + % size of next file in chain curr_block + InitSizeArray(curr_block) = min(nruns(curr_block)-draw_iter,MAX_nruns); + % initialization of next file if necessary + if InitSizeArray(curr_block) + x2 = zeros(InitSizeArray(curr_block),npar); + logpo2 = zeros(InitSizeArray(curr_block),1); + NewFile(curr_block) = NewFile(curr_block) + 1; + draw_index_current_file = 0; + end + end + draw_iter=draw_iter+1; + draw_index_current_file = draw_index_current_file + 1; + end% End of the simulations for one mh-block. + record.AcceptanceRatio(curr_block) = accepted_draws_this_chain/(draw_iter-1); + record.FunctionEvalPerIteration(curr_block) = feval_this_chain/(draw_iter-1); + dyn_waitbar_close(hh); + [record.LastSeeds(curr_block).Unifor, record.LastSeeds(curr_block).Normal] = get_dynare_random_generator_state(); + OutputFileName(block_iter,:) = {[MetropolisFolder,filesep], [ModelName '_mh*_blck' int2str(curr_block) '.mat']}; +end% End of the loop over the mh-blocks. + + +myoutput.record = record; +myoutput.irun = draw_index_current_file; +myoutput.NewFile = NewFile; +myoutput.OutputFileName = OutputFileName; \ No newline at end of file diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m new file mode 100644 index 000000000..2e3804453 --- /dev/null +++ b/matlab/posterior_sampler_initialization.m @@ -0,0 +1,384 @@ +function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... + posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) +%function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, 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] (nblck*npar) vector of starting points for different chains +% o ilogpo2 [double] (nblck*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 fblck [scalar] number of the first MH chain to be run (not equal to 1 in case of recovery) +% o fline [double] (nblck*1) vector of first draw in each chain (not equal to 1 in case of recovery) +% o npar [scalar] number of parameters estimated +% o nblck [scalar] Number of MCM chains requested +% o nruns [double] (nblck*1) number of draws in each chain +% o NewFile [scalar] (nblck*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 (C) 2006-2015 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 = []; +fblck = []; +fline = []; +npar = []; +nblck = []; +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]; + +nblck = options_.mh_nblck; +nruns = ones(nblck,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 nblck > 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(nblck) '\n']); + fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']); + fprintf(fidlog,' \n'); + % Find initial values for the nblck chains... + if isempty(d), + prior_draw(1); + end + if nblck > 1% Case 1: multiple chains + fprintf(fidlog,[' Initial values of the parameters:\n']); + disp('Estimation::mcmc: Searching for initial values...') + ix2 = zeros(nblck,npar); + ilogpo2 = zeros(nblck,1); + for j=1:nblck + validate = 0; + init_iter = 0; + trial = 1; + while validate == 0 && trial <= 10 + if isempty(d), + candidate = prior_draw(0); + else + candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar); + end + 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)]) + 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.') + return + end + fprintf(fidlog,' \n'); + end + fprintf(fidlog,' \n'); + fblck = 1; + fline = ones(nblck,1); + NewFile = ones(nblck,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.Sampler = options_.posterior_sampling_method; + record.Nblck = nblck; + record.MhDraws = zeros(1,3); + record.MhDraws(1,1) = nruns(1); + record.MhDraws(1,2) = AnticipatedNumberOfFiles; + record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile; + record.AcceptanceRatio = zeros(1,nblck); + record.FunctionEvalPerIteration = NaN(1,nblck); + for j=1:nblck + % 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(nblck,npar); + record.LastLogPost = zeros(nblck,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:nblck, + 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(nblck) '\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 ~= nblck + disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!') + disp(['Estimation::mcmc: You declared ' int2str(nblck) ' 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.' ]) + nblck = past_number_of_blocks; + options_.mh_nblck = nblck; + 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(nblck,1)*LastFileNumber; + fline = ones(nblck,1)*(LastLineNumber+1); + else + NewFile = ones(nblck,1)*(LastFileNumber+1); + fline = ones(nblck,1); + end + ilogpo2 = record.LastLogPost; + ix2 = record.LastParameters; + fblck = 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.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); + nblck = record.Nblck;% Number of "parallel" mcmc chains. + options_.mh_nblck = nblck; + if size(record.MhDraws,1) == 1 + OldMh = 0;% The crashed metropolis was the first session. + else + OldMh = 1;% The crashed metropolis wasn't the first session. + end + % Default initialization: + if OldMh + ilogpo2 = record.LastLogPost; + ix2 = record.LastParameters; + else + ilogpo2 = record.InitialLogPost; + ix2 = record.InitialParameters; + end + % Set NewFile, a nblck*1 vector of integers, and fline (first line), a nblck*1 vector of integers. + if OldMh + 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. + if LastLineNumberInThePreviousMh < MAX_nruns% Test if the last mh files of the previous session are not complete (yes) + NewFile = ones(nblck,1)*LastFileNumberInThePreviousMh; + fline = ones(nblck,1)*(LastLineNumberInThePreviousMh+1); + else% The last mh files of the previous session are complete. + NewFile = ones(nblck,1)*(LastFileNumberInThePreviousMh+1); + fline = ones(nblck,1); + end + else + LastLineNumberInThePreviousMh = 0; + LastFileNumberInThePreviousMh = 0; + NewFile = ones(nblck,1); + fline = ones(nblck,1); + end + % Set fblck (First block), an integer targeting the crashed mcmc chain. + fblck = 1; + % How many mh files should we have ? + ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1); + ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*nblck; + % I count the total number of saved mh files... + AllMhFiles = dir([BaseName '_mh*_blck*.mat']); + TotalNumberOfMhFiles = length(AllMhFiles); + % And I quit if I can't find a crashed mcmc chain. + 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 + % I count the number of saved mh files per block. + NumberOfMhFilesPerBlock = zeros(nblck,1); + for b = 1:nblck + BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']); + NumberOfMhFilesPerBlock(b) = length(BlckMhFiles); + end + % Is there a chain with less mh files than expected ? + while fblck <= nblck + if NumberOfMhFilesPerBlock(fblck) < ExpectedNumberOfMhFilesPerBlock + disp(['Estimation::mcmc: Chain ' int2str(fblck) ' is not complete!']) + break + % The mh_recover session will start from chain fblck. + else + disp(['Estimation::mcmc: Chain ' int2str(fblck) ' is complete!']) + end + fblck = fblck+1; + end + % How many mh-files are saved in this block? + NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(fblck); + % 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 OldMh + NumberOfSavedMhFilesInTheCrashedBlck = NumberOfSavedMhFilesInTheCrashedBlck - LastFileNumberInThePreviousMh; + end + NumberOfSavedMhFiles = NumberOfSavedMhFilesInTheCrashedBlck+LastFileNumberInThePreviousMh; + % Correct initial conditions. + if NumberOfSavedMhFiles + load([BaseName '_mh' int2str(NumberOfSavedMhFiles) '_blck' int2str(fblck) '.mat']); + ilogpo2(fblck) = logpo2(end); + ix2(fblck,:) = x2(end,:); + end +end \ No newline at end of file diff --git a/matlab/posterior_sampler_iteration.m b/matlab/posterior_sampler_iteration.m new file mode 100644 index 000000000..bbd642d75 --- /dev/null +++ b/matlab/posterior_sampler_iteration.m @@ -0,0 +1,93 @@ +function [par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun,last_draw, last_posterior, sampler_options,varargin) + +% function [par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun,last_draw, last_posterior, sampler_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_) +% 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 + +% Copyright (C) 2015 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 . + + +posterior_sampling_method = sampler_options.posterior_sampling_method; +mh_bounds = sampler_options.bounds; + +switch posterior_sampling_method + case 'slice' + + [par, logpost, neval] = slice_sampler(TargetFun,last_draw, [mh_bounds.lb mh_bounds.ub], sampler_options,varargin{:}); + accepted = 1; + case 'random_walk_metropolis_hastings' + neval = 1; + ProposalFun = sampler_options.proposal_distribution; + proposal_covariance_Cholesky_decomposition = sampler_options.proposal_covariance_Cholesky_decomposition; + n = sampler_options.n; + + par = feval(ProposalFun, last_draw, proposal_covariance_Cholesky_decomposition, n); + if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub ) + try + logpost = - feval(TargetFun, par(:),varargin{:}); + catch + logpost = -inf; + end + else + logpost = -inf; + end + r = logpost-last_posterior; + if (logpost > -inf) && (log(rand) < r) + accepted = 1; + else + accepted = 0; + par = last_draw; + logpost = last_posterior; + end + case 'independent_metropolis_hastings' + neval = 1; + ProposalFun = sampler_options.proposal_distribution; + ProposalDensity = sampler_options.ProposalDensity; + proposal_covariance_Cholesky_decomposition = sampler_options.proposal_covariance_Cholesky_decomposition; + n = sampler_options.n; + xparam1 = sampler_options.xparam1; + par = feval(ProposalFun, sampler_options.xparam1, proposal_covariance_Cholesky_decomposition, n); + if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub ) + try + logpost = - feval(TargetFun, par(:),varargin{:}); + catch + logpost = -inf; + end + else + logpost = -inf; + end + r = logpost - last_posterior + ... + log(feval(ProposalDensity, last_draw, xparam1, proposal_covariance_Cholesky_decomposition, n)) - ... + log(feval(ProposalDensity, par, xparam1, proposal_covariance_Cholesky_decomposition, n)); + if (logpost > -inf) && (log(rand) < r) + accepted = 1; + else + accepted = 0; + par = last_draw; + logpost = last_posterior; + end +end \ No newline at end of file