1function SampleAddress = selec_posterior_draws(SampleSize,drsize) 2% Selects a sample of draws from the posterior distribution and if nargin>1 3% saves the draws in _pdraws mat files (metropolis folder). If drsize>0 4% the dr structure, associated to the parameters, is also saved in _pdraws. 5% This routine is more efficient than metropolis_draw.m because here an 6% _mh file cannot be opened twice. 7% 8% INPUTS 9% o SampleSize [integer] Size of the sample to build. 10% o drsize [double] structure dr is drsize megaoctets. 11% 12% OUTPUTS 13% o SampleAddress [integer] A (SampleSize*4) array, each line specifies the 14% location of a posterior draw: 15% Column 2 --> Chain number 16% Column 3 --> (mh) File number 17% Column 4 --> (mh) line number 18% 19% SPECIAL REQUIREMENTS 20% None. 21% 22 23% Copyright (C) 2006-2017 Dynare Team 24% 25% This file is part of Dynare. 26% 27% Dynare is free software: you can redistribute it and/or modify 28% it under the terms of the GNU General Public License as published by 29% the Free Software Foundation, either version 3 of the License, or 30% (at your option) any later version. 31% 32% Dynare is distributed in the hope that it will be useful, 33% but WITHOUT ANY WARRANTY; without even the implied warranty of 34% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 35% GNU General Public License for more details. 36% 37% You should have received a copy of the GNU General Public License 38% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 39 40global M_ options_ estim_params_ oo_ 41 42% Number of parameters: 43npar = estim_params_.nvx; 44npar = npar + estim_params_.nvn; 45npar = npar + estim_params_.ncx; 46npar = npar + estim_params_.ncn; 47npar = npar + estim_params_.np; 48 49% Select one task: 50switch nargin 51 case 1 52 info = 0; 53 case 2 54 MAX_mega_bytes = 10;% Should be an option... 55 if drsize>0 56 info=2; 57 else 58 info=1; 59 end 60 drawsize = drsize+npar*8/1048576; 61 otherwise 62 error(['selec_posterior_draws:: Unexpected number of input arguments!']) 63end 64 65MetropolisFolder = CheckPath('metropolis',M_.dname); 66ModelName = M_.fname; 67BaseName = [MetropolisFolder filesep ModelName]; 68 69% Get informations about the mcmc: 70load_last_mh_history_file(MetropolisFolder, ModelName); 71FirstMhFile = record.KeepedDraws.FirstMhFile; 72FirstLine = record.KeepedDraws.FirstLine; 73TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); 74LastMhFile = TotalNumberOfMhFiles; 75TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); 76NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); 77MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); 78mh_nblck = options_.mh_nblck; 79 80% Randomly select draws in the posterior distribution: 81SampleAddress = zeros(SampleSize,4); 82for i = 1:SampleSize 83 ChainNumber = ceil(rand*mh_nblck); 84 DrawNumber = ceil(rand*NumberOfDraws); 85 SampleAddress(i,1) = DrawNumber; 86 SampleAddress(i,2) = ChainNumber; 87 if DrawNumber <= MAX_nruns-FirstLine+1 88 MhFileNumber = FirstMhFile; 89 MhLineNumber = FirstLine+DrawNumber-1; 90 else 91 DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1); 92 MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns); 93 MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns; 94 end 95 SampleAddress(i,3) = MhFileNumber; 96 SampleAddress(i,4) = MhLineNumber; 97end 98SampleAddress = sortrows(SampleAddress,[3 2]); 99 100% Selected draws in the posterior distribution, and if drsize>0 101% reduced form solutions, are saved on disk. 102if info 103 %delete old stale files before creating new ones 104 delete_stale_file([BaseName '_posterior_draws*.mat']) 105 if SampleSize*drawsize <= MAX_mega_bytes% The posterior draws are saved in one file. 106 pdraws = cell(SampleSize,info); 107 old_mhfile = 0; 108 old_mhblck = 0; 109 for i = 1:SampleSize 110 mhfile = SampleAddress(i,3); 111 mhblck = SampleAddress(i,2); 112 if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck) 113 load([BaseName '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2') 114 end 115 pdraws(i,1) = {x2(SampleAddress(i,4),:)}; 116 if info-1 117 set_parameters(pdraws{i,1}); 118 [dr,info,M_,options_,oo_] =compute_decision_rules(M_,options_,oo_); 119 pdraws(i,2) = { dr }; 120 end 121 old_mhfile = mhfile; 122 old_mhblck = mhblck; 123 end 124 clear('x2') 125 save([BaseName '_posterior_draws1.mat'],'pdraws','estim_params_') 126 else% The posterior draws are saved in xx files. 127 NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize); 128 NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes); 129 NumberOfLines = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile; 130 linee = 0; 131 fnum = 1; 132 pdraws = cell(NumberOfDrawsPerFile,info); 133 old_mhfile = 0; 134 old_mhblck = 0; 135 for i=1:SampleSize 136 linee = linee+1; 137 mhfile = SampleAddress(i,3); 138 mhblck = SampleAddress(i,2); 139 if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck) 140 load([BaseName '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2') 141 end 142 pdraws(linee,1) = {x2(SampleAddress(i,4),:)}; 143 if info-1 144 set_parameters(pdraws{linee,1}); 145 [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_); 146 pdraws(linee,2) = { dr }; 147 end 148 old_mhfile = mhfile; 149 old_mhblck = mhblck; 150 if fnum < NumberOfFiles && linee == NumberOfDrawsPerFile 151 linee = 0; 152 save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_') 153 fnum = fnum+1; 154 if fnum < NumberOfFiles 155 pdraws = cell(NumberOfDrawsPerFile,info); 156 else 157 pdraws = cell(NumberOfLines,info); 158 end 159 end 160 end 161 save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_') 162 end 163end