From a6281ac9f95e06bbf5819710eaeeba7d92b1aa04 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Sun, 24 Jul 2011 19:50:56 +0200 Subject: [PATCH] adding a first attempt for adaptive Metropolis Hastings --- matlab/adaptive_metropolis_hastings.m | 132 +++++++++++++++++++++++ matlab/random_walk_metropolis_hastings.m | 6 +- 2 files changed, 135 insertions(+), 3 deletions(-) create mode 100644 matlab/adaptive_metropolis_hastings.m diff --git a/matlab/adaptive_metropolis_hastings.m b/matlab/adaptive_metropolis_hastings.m new file mode 100644 index 000000000..1875dac9e --- /dev/null +++ b/matlab/adaptive_metropolis_hastings.m @@ -0,0 +1,132 @@ +function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin) +%function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin) +% Random walk Metropolis-Hastings algorithm. +% +% 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 varargin list of argument following mh_bounds +% +% OUTPUTS +% None +% +% ALGORITHM +% Metropolis-Hastings. +% +% SPECIAL REQUIREMENTS +% None. +% +% PARALLEL CONTEXT +% The most computationally intensive part of this function may be executed +% in parallel. The code sutable to be executed in +% parallel on multi core or cluster machine (in general a 'for' cycle), +% is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion. +% Then the DYNARE parallel package contain a set of pairs matlab functions that can be executed in +% parallel and called name_function.m and name_function_core.m. +% In addition in parallel package we have second set of functions used +% to manage the parallel computation. +% +% 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 funtions. + +% Copyright (C) 2006-2011 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 . + +global M_ options_ bayestopt_ estim_params_ oo_ + + +old_options = options_; + +accept_target = options_.amh.accept_target; + +options_.mh_jscale = tune_scale_parameter(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin{:}); + +for i=1:options_.amh.cova_steps + options_.mh_replic = options_.amh.cova_replic; + random_walk_metropolis_hastings(TargetFun,ProposalFun, ... + xparam1,vv,mh_bounds,varargin{:}); + CutSample(M_,options_,estim_params_); + [junk,vv] = compute_mh_covariance_matrix(); + bayestopt_.jscale = tune_scale_parameter(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin{:}); +end + +options_.mh_replic = old_options.mh_replic; +record = random_walk_metropolis_hastings(TargetFun,ProposalFun, ... + xparam1,vv,mh_bounds,varargin{:}); + + + +function selected_scale = tune_scale_parameter(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin) +global options_ bayestopt_ + +selected_scale = []; + +maxit = options_.amh.scale_tuning_maxit; +accept_target = options_.amh.accept_target; +test_runs = options_.amh.scale_tuning_test_runs; +tolerance = options_.amh.scale_tuning_tolerance; +Scales = zeros(maxit,1); +AvRates = zeros(maxit,1); +Scales(1) = options_.mh_jscale; + +for i=1:maxit + options_.mh_replic = options_.amh.scale_tuning_blocksize; + bayestopt_.jscale = Scales(i); + record = random_walk_metropolis_hastings(TargetFun,ProposalFun, ... + xparam1,vv, ... + mh_bounds,varargin{:}); + AvRates(i) = mean(record.AcceptationRates); + disp(AvRates(1:i)') + + if i >= test_runs + a_mean = mean(AvRates((i-test_runs+1):i)); + if (a_mean > (1-tolerance)*accept_target) && ... + (a_mean < (1+tolerance)*accept_target) + selected_scale = mean(Scales((i-test_runs+1):i)); + disp(['Selected scale: ' num2str(selected_scale)]) + return + end + end + + if i == 1 + if AvRates(1) > accept_target + Scales(i+1) = 2*Scales(i); + else + Scales(i+1) = Scales(i)/2; + end + elseif i < maxit + X = [ones(i,1) Scales(1:i)]; + b = X\(AvRates(1:i)-accept_target); + Scales(i+1) = -b(1)/b(2); + if Scales(i+1) < 0.05 + Scales(i+1) = 0.05; + elseif Scales(i+1) > 2 + Scales(i+1) = 2; + end + else + error('AMH scale tuning: tuning didn''t converge') + end + + options_.load_mh_file = 1; + + disp(Scales(1:i)') +end \ No newline at end of file diff --git a/matlab/random_walk_metropolis_hastings.m b/matlab/random_walk_metropolis_hastings.m index cf55d0e65..48c20fa94 100644 --- a/matlab/random_walk_metropolis_hastings.m +++ b/matlab/random_walk_metropolis_hastings.m @@ -1,5 +1,5 @@ -function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin) -%function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin) +function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin) +%function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin) % Random walk Metropolis-Hastings algorithm. % % INPUTS @@ -11,7 +11,7 @@ function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou % o varargin list of argument following mh_bounds % % OUTPUTS -% None +% o record [struct] structure describing the iterations % % ALGORITHM % Metropolis-Hastings.