metropolis_draw

 function [xparams, logpost]=metropolis_draw(init) 
 Builds draws from metropolis

 INPUTS:
   init:              scalar equal to 1 (first call) or 0

 OUTPUTS:
   xparams:           vector of estimated parameters
   logpost:           log of posterior density
   
 SPECIAL REQUIREMENTS
   none

0001 function [xparams, logpost]=metropolis_draw(init)
0002 % function [xparams, logpost]=metropolis_draw(init)
0003 % Builds draws from metropolis
0004 %
0005 % INPUTS:
0006 %   init:              scalar equal to 1 (first call) or 0
0007 %
0008 % OUTPUTS:
0009 %   xparams:           vector of estimated parameters
0010 %   logpost:           log of posterior density
0011 %
0012 % SPECIAL REQUIREMENTS
0013 %   none
0014 
0015 % Copyright (C) 2003-2009 Dynare Team
0016 %
0017 % This file is part of Dynare.
0018 %
0019 % Dynare is free software: you can redistribute it and/or modify
0020 % it under the terms of the GNU General Public License as published by
0021 % the Free Software Foundation, either version 3 of the License, or
0022 % (at your option) any later version.
0023 %
0024 % Dynare is distributed in the hope that it will be useful,
0025 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0026 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0027 % GNU General Public License for more details.
0028 %
0029 % You should have received a copy of the GNU General Public License
0030 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0031 
0032 global options_ estim_params_ M_
0033 persistent mh_nblck NumberOfDraws fname FirstLine FirstMhFile MAX_nruns
0034 
0035 if init
0036     nvx  = estim_params_.nvx;
0037     nvn  = estim_params_.nvn;
0038     ncx  = estim_params_.ncx;
0039     ncn  = estim_params_.ncn;
0040     np   = estim_params_.np ;
0041     npar = nvx+nvn+ncx+ncn+np;
0042     MhDirectoryName = CheckPath('metropolis',M_.dname);
0043     fname = [ MhDirectoryName '/' M_.fname];
0044     load([ fname '_mh_history']);
0045     FirstMhFile = record.KeepedDraws.FirstMhFile;
0046     FirstLine = record.KeepedDraws.FirstLine; 
0047     TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); 
0048     LastMhFile = TotalNumberOfMhFiles; 
0049     TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
0050     NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
0051     MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
0052     mh_nblck = options_.mh_nblck;
0053     return
0054 end
0055 
0056 ChainNumber = ceil(rand*mh_nblck);
0057 DrawNumber  = ceil(rand*NumberOfDraws);
0058 
0059 if DrawNumber <= MAX_nruns-FirstLine+1
0060     MhFilNumber = FirstMhFile;
0061     MhLine = FirstLine+DrawNumber-1;
0062 else
0063     DrawNumber  = DrawNumber-(MAX_nruns-FirstLine+1);
0064     MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns); 
0065     MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns;
0066 end
0067 
0068 load( [ fname '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2');
0069 xparams = x2(MhLine,:);
0070 logpost= logpo2(MhLine);