2008-06-17 23:23:02 +02:00
|
|
|
function SampleAddress = selec_posterior_draws(SampleSize,drsize)
|
|
|
|
% Selects a sample of draws from the posterior distribution and if nargin>1
|
|
|
|
% saves the draws in _pdraws mat files (metropolis folder). If drsize>0
|
2011-09-17 12:33:07 +02:00
|
|
|
% the dr structure, associated to the parameters, is also saved in _pdraws.
|
2008-06-17 23:23:02 +02:00
|
|
|
% This routine is more efficient than metropolis_draw.m because here an
|
2011-09-17 12:33:07 +02:00
|
|
|
% _mh file cannot be opened twice.
|
|
|
|
%
|
2006-10-10 21:37:26 +02:00
|
|
|
% INPUTS
|
2008-06-17 23:23:02 +02:00
|
|
|
% o SampleSize [integer] Size of the sample to build.
|
2011-09-17 12:33:07 +02:00
|
|
|
% o drsize [double] structure dr is drsize megaoctets.
|
2008-01-17 16:52:13 +01:00
|
|
|
%
|
2007-06-16 10:59:31 +02:00
|
|
|
% OUTPUTS
|
2011-09-17 12:33:07 +02:00
|
|
|
% o SampleAddress [integer] A (SampleSize*4) array, each line specifies the
|
|
|
|
% location of a posterior draw:
|
2006-10-10 21:37:26 +02:00
|
|
|
% Column 2 --> Chain number
|
2011-09-17 12:33:07 +02:00
|
|
|
% Column 3 --> (mh) File number
|
2006-10-10 21:37:26 +02:00
|
|
|
% Column 4 --> (mh) line number
|
|
|
|
%
|
|
|
|
% SPECIAL REQUIREMENTS
|
|
|
|
% None.
|
2011-09-17 12:33:07 +02:00
|
|
|
%
|
2008-08-01 14:40:33 +02:00
|
|
|
|
2011-02-10 17:23:22 +01:00
|
|
|
% Copyright (C) 2006-2011 Dynare Team
|
2008-08-01 14:40:33 +02:00
|
|
|
%
|
|
|
|
% 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 <http://www.gnu.org/licenses/>.
|
2008-06-17 23:23:02 +02:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
global M_ options_ estim_params_ oo_
|
|
|
|
|
|
|
|
% Number of parameters:
|
|
|
|
npar = estim_params_.nvx;
|
|
|
|
npar = npar + estim_params_.nvn;
|
|
|
|
npar = npar + estim_params_.ncx;
|
|
|
|
npar = npar + estim_params_.ncn;
|
|
|
|
npar = npar + estim_params_.np;
|
|
|
|
|
2011-09-17 12:33:07 +02:00
|
|
|
% Select one task:
|
2009-12-16 18:17:34 +01:00
|
|
|
switch nargin
|
|
|
|
case 1
|
|
|
|
info = 0;
|
|
|
|
case 2
|
2010-08-19 15:39:33 +02:00
|
|
|
MAX_mega_bytes = 10;% Should be an option...
|
2009-12-16 18:17:34 +01:00
|
|
|
if drsize>0
|
|
|
|
info=2;
|
|
|
|
else
|
|
|
|
info=1;
|
2008-06-17 23:23:02 +02:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
drawsize = drsize+npar*8/1048576;
|
|
|
|
otherwise
|
|
|
|
error(['selec_posterior_draws:: Unexpected number of input arguments!'])
|
|
|
|
end
|
|
|
|
|
2011-09-17 12:33:07 +02:00
|
|
|
% Get informations about the mcmc:
|
2011-12-15 11:56:23 +01:00
|
|
|
MhDirectoryName = CheckPath('metropolis',M_.dname);
|
2009-12-16 18:17:34 +01:00
|
|
|
fname = [ MhDirectoryName '/' M_.fname];
|
|
|
|
load([ fname '_mh_history.mat']);
|
|
|
|
FirstMhFile = record.KeepedDraws.FirstMhFile;
|
2011-09-17 12:33:07 +02:00
|
|
|
FirstLine = record.KeepedDraws.FirstLine;
|
|
|
|
TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
|
|
|
|
LastMhFile = TotalNumberOfMhFiles;
|
2009-12-16 18:17:34 +01:00
|
|
|
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
|
|
|
|
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
|
|
|
|
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
|
|
|
|
mh_nblck = options_.mh_nblck;
|
|
|
|
|
|
|
|
% Randomly select draws in the posterior distribution:
|
|
|
|
SampleAddress = zeros(SampleSize,4);
|
|
|
|
for i = 1:SampleSize
|
|
|
|
ChainNumber = ceil(rand*mh_nblck);
|
|
|
|
DrawNumber = ceil(rand*NumberOfDraws);
|
|
|
|
SampleAddress(i,1) = DrawNumber;
|
|
|
|
SampleAddress(i,2) = ChainNumber;
|
|
|
|
if DrawNumber <= MAX_nruns-FirstLine+1
|
|
|
|
MhFileNumber = FirstMhFile;
|
|
|
|
MhLineNumber = FirstLine+DrawNumber-1;
|
|
|
|
else
|
|
|
|
DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1);
|
2011-09-17 12:33:07 +02:00
|
|
|
MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
|
2009-12-16 18:17:34 +01:00
|
|
|
MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns;
|
2007-01-22 17:34:54 +01:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
SampleAddress(i,3) = MhFileNumber;
|
|
|
|
SampleAddress(i,4) = MhLineNumber;
|
|
|
|
end
|
|
|
|
SampleAddress = sortrows(SampleAddress,[3 2]);
|
|
|
|
|
|
|
|
% Selected draws in the posterior distribution, and if drsize>0
|
|
|
|
% reduced form solutions, are saved on disk.
|
|
|
|
if info
|
|
|
|
if SampleSize*drawsize <= MAX_mega_bytes% The posterior draws are saved in one file.
|
|
|
|
pdraws = cell(SampleSize,info);
|
|
|
|
old_mhfile = 0;
|
|
|
|
old_mhblck = 0;
|
|
|
|
for i = 1:SampleSize
|
|
|
|
mhfile = SampleAddress(i,3);
|
|
|
|
mhblck = SampleAddress(i,2);
|
2011-02-10 17:23:22 +01:00
|
|
|
if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
|
2009-12-16 18:17:34 +01:00
|
|
|
load([fname '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
|
2008-06-17 23:23:02 +02:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
pdraws(i,1) = {x2(SampleAddress(i,4),:)};
|
|
|
|
if info-1
|
|
|
|
set_parameters(pdraws{i,1});
|
2011-09-17 12:33:07 +02:00
|
|
|
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
|
2009-12-16 18:17:34 +01:00
|
|
|
pdraws(i,2) = { dr };
|
|
|
|
end
|
|
|
|
old_mhfile = mhfile;
|
|
|
|
old_mhblck = mhblck;
|
|
|
|
end
|
|
|
|
clear('x2')
|
2011-04-06 13:08:22 +02:00
|
|
|
save([fname '_posterior_draws1.mat'],'pdraws')
|
2009-12-16 18:17:34 +01:00
|
|
|
else% The posterior draws are saved in xx files.
|
|
|
|
NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize);
|
2011-09-17 12:33:07 +02:00
|
|
|
NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes);
|
2009-12-16 18:17:34 +01:00
|
|
|
NumberOfLines = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile;
|
|
|
|
linee = 0;
|
|
|
|
fnum = 1;
|
|
|
|
pdraws = cell(NumberOfDrawsPerFile,info);
|
|
|
|
old_mhfile = 0;
|
|
|
|
old_mhblck = 0;
|
|
|
|
for i=1:SampleSize
|
|
|
|
linee = linee+1;
|
|
|
|
mhfile = SampleAddress(i,3);
|
|
|
|
mhblck = SampleAddress(i,2);
|
2011-02-10 17:23:22 +01:00
|
|
|
if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
|
2009-12-16 18:17:34 +01:00
|
|
|
load([fname '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
|
|
|
|
end
|
|
|
|
pdraws(linee,1) = {x2(SampleAddress(i,4),:)};
|
|
|
|
if info-1
|
|
|
|
set_parameters(pdraws{linee,1});
|
2011-09-17 12:33:07 +02:00
|
|
|
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
|
2009-12-16 18:17:34 +01:00
|
|
|
pdraws(linee,2) = { dr };
|
|
|
|
end
|
|
|
|
old_mhfile = mhfile;
|
|
|
|
old_mhblck = mhblck;
|
|
|
|
if fnum < NumberOfFiles && linee == NumberOfDrawsPerFile
|
|
|
|
linee = 0;
|
2011-04-06 13:08:22 +02:00
|
|
|
save([fname '_posterior_draws' num2str(fnum) '.mat'],'pdraws')
|
2009-12-16 18:17:34 +01:00
|
|
|
fnum = fnum+1;
|
|
|
|
if fnum < NumberOfFiles
|
|
|
|
pdraws = cell(NumberOfDrawsPerFile,info);
|
|
|
|
else
|
|
|
|
pdraws = cell(NumberOfLines,info);
|
2008-06-17 23:23:02 +02:00
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
2011-04-06 13:08:22 +02:00
|
|
|
save([fname '_posterior_draws' num2str(fnum) '.mat'],'pdraws')
|
2009-12-16 18:17:34 +01:00
|
|
|
end
|
|
|
|
end
|