Commit 433a6816 authored by Johannes Pfeifer 's avatar Johannes Pfeifer Committed by Stéphane Adjemian

TaRB: fix bug where incorrect last posterior was returned if last draw was rejected

parent 157641ac
...@@ -4,16 +4,29 @@ function [par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFu ...@@ -4,16 +4,29 @@ function [par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFu
% posterior samplers % posterior samplers
% %
% INPUTS % INPUTS
% posterior_sampler_options: posterior sampler options % TargetFun: string storing the objective function (e.g. 'dsge_likelihood.m')
% options_: structure storing the options % last_draw: parameter vector in last iteration
% last_posterior: value of the posterior in last iteration
% sampler_options: posterior sampler options
% dataset_: the dataset after required transformation
% dataset_info: Various informations about the dataset (descriptive statistics and missing observations).
% options_: structure storing the options
% M_: structure storing the model information
% estim_params_: structure storing information about estimated parameters
% bayestopt_: structure storing information about priors
% mh_bounds: structure containing prior bounds
% oo_: structure storing the results
%
% OUTPUTS % OUTPUTS
% posterior_sampler_options: checked posterior sampler options % par: last accepted parameter vector
% logpost: value of the posterior after current iteration
% accepted: share of proposed draws that were accepted
% neval: number of evaluations (>1 only for slice)
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright (C) 2015-16 Dynare Team % Copyright (C) 2015-18 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -121,28 +134,31 @@ switch posterior_sampling_method ...@@ -121,28 +134,31 @@ switch posterior_sampling_method
else else
logpost = -inf; logpost = -inf;
end end
%get ratio of proposal densities, required because proposal depends
%on current mode via Hessian and is thus not symmetric anymore if (logpost > -inf)
if strcmpi(sampler_options.proposal_distribution,'rand_multivariate_normal') %get ratio of proposal densities, required because proposal depends
proposal_density_proposed_move_forward=multivariate_normal_pdf(proposed_par,xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n); %on current mode via Hessian and is thus not symmetric anymore
proposal_density_proposed_move_backward=multivariate_normal_pdf(par_start_current_block',xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n); if strcmpi(sampler_options.proposal_distribution,'rand_multivariate_normal')
elseif strcmpi(sampler_options.proposal_distribution,'rand_multivariate_student') proposal_density_proposed_move_forward=multivariate_normal_pdf(proposed_par,xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n);
proposal_density_proposed_move_forward=multivariate_student_pdf(proposed_par,xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n); proposal_density_proposed_move_backward=multivariate_normal_pdf(par_start_current_block',xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n);
proposal_density_proposed_move_backward=multivariate_student_pdf(par_start_current_block',xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n); elseif strcmpi(sampler_options.proposal_distribution,'rand_multivariate_student')
end proposal_density_proposed_move_forward=multivariate_student_pdf(proposed_par,xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n);
accprob=logpost-last_posterior+ log(proposal_density_proposed_move_backward)-log(proposal_density_proposed_move_forward); %Formula (6), Chib/Ramamurthy proposal_density_proposed_move_backward=multivariate_student_pdf(par_start_current_block',xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n);
end
if (logpost > -inf) && (log(rand) < accprob) accprob=logpost-last_posterior+ log(proposal_density_proposed_move_backward)-log(proposal_density_proposed_move_forward); %Formula (6), Chib/Ramamurthy
current_draw(indices(blocks==block_iter,1))=proposed_par; if (log(rand) < accprob)
last_posterior=logpost; current_draw(indices(blocks==block_iter,1))=proposed_par;
accepted_draws_counter =accepted_draws_counter +1; last_posterior=logpost;
else %no updating accepted_draws_counter =accepted_draws_counter +1;
%do nothing, keep old value else %no updating
%do nothing, keep old value
end
end end
end end
accepted=accepted_draws_counter/blocked_draws_counter; accepted=accepted_draws_counter/blocked_draws_counter;
par = current_draw; par = current_draw;
neval=1; neval=1;
logpost = last_posterior; %make sure not a temporary draw is returned;
case 'independent_metropolis_hastings' case 'independent_metropolis_hastings'
neval = 1; neval = 1;
ProposalFun = sampler_options.proposal_distribution; ProposalFun = sampler_options.proposal_distribution;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment