From f2f072bce13525abb3b07062943b209b3ed25ab5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Scylla=29?= Date: Tue, 26 Nov 2013 11:28:40 +0100 Subject: [PATCH] Added options to the internals command to load mh-history files or display informations about MCMC settings and state. For instance, assuming that fs2000a.mod is in the current Matlab/Octave's directory and that a metropolis has been run, the following command: >> internals --display-mh-history fs2000a produces the following output in the Matlab/Octave's command window: ================================ MCMC set-up for fs2000a mod file ================================ MCMC chain number 1: -------------------- o Number of MCMC files is 1 o Number of draws is 6000 o Acceptance ratio is 35.28% o Last value of the posterior kernel is: 1221.81656 o State of the chain: || Initial | Current ++++++++++++++++++++++++ e_a || 0.01570 | 0.01608 e_m || 0.00528 | 0.00512 alp || 0.35026 | 0.36909 bet || 0.99300 | 0.99308 gam || 0.00086 | 0.00038 mst || 1.00088 | 1.00147 rho || 0.65926 | 0.71132 psi || 0.66428 | 0.60611 del || 0.00977 | 0.00198 MCMC chain number 2: -------------------- o Number of MCMC files is 1 o Number of draws is 6000 o Acceptance ratio is 33.38% o Last value of the posterior kernel is: 1220.39229 o State of the chain: || Initial | Current +++++++++++++++++++++++++ e_a || 0.01664 | 0.01576 e_m || 0.00489 | 0.00504 alp || 0.39701 | 0.36630 bet || 0.99305 | 0.99219 gam || 0.00067 | -0.00025 mst || 1.00186 | 1.00218 rho || 0.70060 | 0.63847 psi || 0.65623 | 0.74668 del || 0.00666 | 0.00842 while the command: >> internals --load-mh-history loads the content of the record structure (saved in the last mh-history file) in Matlab/Octave's workspace under the name mcmc_informations. --- matlab/internals.m | 81 +++++++++++++++++++++++++++++ matlab/load_first_mh_history_file.m | 39 ++++++++++++++ 2 files changed, 120 insertions(+) create mode 100644 matlab/load_first_mh_history_file.m diff --git a/matlab/internals.m b/matlab/internals.m index 98225669d..02205b4f4 100644 --- a/matlab/internals.m +++ b/matlab/internals.m @@ -90,6 +90,87 @@ if strcmpi(flag,'--test') return end +if strcmpi(flag,'--load-mh-history') || strcmpi(flag,'--display-mh-history') + switch length(varargin) + case 3 + fname = varargin{1}; + if ~isequal(varargin{2},'in') + error('internals:: Calling sequence must be of the form: internals --load-mh-history fname in dname') + end + dname = varargin{3}; + case 1 + fname = varargin{1}; + dname = varargin{1}; + otherwise + error('internals:: Wrong calling sequence! You should read the manual...') + end + o = load_last_mh_history_file([dname filesep 'metropolis'],fname); + if strcmpi(flag,'--load-mh-history') + assignin('caller','mcmc_informations',o.record); + else + oo = load_first_mh_history_file([dname filesep 'metropolis'],fname); + local = load([fname '_results'],'bayestopt_'); + names = local.bayestopt_.name; %evalin('base','bayestopt_.name'); + str = ['MCMC set-up for ' fname ' mod file']; + ltr = length(str); + skipline() + disp(repmat('=',1,ltr)) + disp(str) + disp(repmat('=',1,ltr)) + skipline(2) + for b=1:o.Nblck + str = ['MCMC chain number ' num2str(b) ':']; + ltr = length(str); + disp(str); + disp(repmat('-',1,ltr)); + skipline() + disp([' o Number of MCMC files is ' num2str(sum(o.MhDraws(:,2)))]); + disp([' o Number of draws is ' num2str(sum(o.MhDraws(:,1)))]); + disp([' o Acceptance ratio is ' num2str(o.AcceptationRates(b)*100,'%6.2f') '%']); + disp([' o Last value of the posterior kernel is: ' num2str(o.LastLogPost(b),'%10.5f')]) + disp([' o State of the chain:']) + skipline() + d1 = num2str(transpose(oo.InitialParameters(b,:)),'%10.5f\n'); + d2 = num2str(transpose(o.LastParameters(b,:)),'%10.5f\n'); + d1s = size(d1,2); + d2s = size(d2,2); + c0 = repmat(' ',length(names)+2,1); + c1 = char(' ', repmat('+',1,size(cell2mat(names),2)), cell2mat(names)); + s1 = char(' || ','++++',repmat(' || ', length(names),1)); + t1 = repmat(' ',1,d1s); + if d1s<=7 + t1 = 'Initial'; + else + diff = d1s-7; + if isequal(mod(diff,2),0) + start = diff/2+1; + else + start = (diff-1)/2+1; + end + t1(start:start+6) = 'Initial'; + end + c2 = char(t1,repmat('+',1,size(d1,2)),d1); + s2 = char(' | ','+++',repmat(' | ', length(names),1)); + t2 = repmat(' ',1,d2s); + if d2s<=7 + t2 = 'Current'; + else + diff = d2s-7; + if isequal(mod(diff,2),0) + start = diff/2+1; + else + start = (diff-1)/2+1; + end + t2(start:start+6) = 'Current'; + end + c3 = char(t2,repmat('+',1,size(d2,2)), d2); + disp([c0, c1, s1, c2, s2, c3]); + skipline() + end + end + return +end + if strcmpi(flag,'--info') if nargin==2 dynare_config([],0); diff --git a/matlab/load_first_mh_history_file.m b/matlab/load_first_mh_history_file.m new file mode 100644 index 000000000..3177fafe4 --- /dev/null +++ b/matlab/load_first_mh_history_file.m @@ -0,0 +1,39 @@ +function info = load_first_mh_history_file(MetropolisFolder, ModelName) + +% This routine requires that the MCMC draws were obtained with a dynare version greater than 4.3.3. + +% Copyright (C) 2013 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 . + +% record is also a Matlab function. +record = 0; + +% Get the list of all the mh_history files. +BaseName = [MetropolisFolder filesep ModelName]; +mh_history_files = dir([BaseName '_mh_history_*.mat']); + +if isequal(length(mh_history_files),0) + error(['Estimation::load_mh_file: I cannot find any mh-history file in ' MetropolisFolder '!']) +end + +load([BaseName '_mh_history_0.mat']); + +if isequal(nargout,0) + assignin('caller', 'record', record); +else + info = record; +end \ No newline at end of file