selec_posterior_draws

 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
 the dr structure, associated to the parameters, is also saved in _pdraws.
 This routine is more efficient than metropolis_draw.m because here an
 _mh file cannot be opened twice.

 INPUTS
   o SampleSize     [integer]  Size of the sample to build.
   o drsize         [double]   structure dr is drsize megaoctets.

 OUTPUTS
   o SampleAddress  [integer]  A (SampleSize*4) array, each line specifies the
                               location of a posterior draw:
                                  Column 2 --> Chain number
                                  Column 3 --> (mh) File number
                                  Column 4 --> (mh) line number

 SPECIAL REQUIREMENTS
   None.

0001 function SampleAddress = selec_posterior_draws(SampleSize,drsize)
0002 % Selects a sample of draws from the posterior distribution and if nargin>1
0003 % saves the draws in _pdraws mat files (metropolis folder). If drsize>0
0004 % the dr structure, associated to the parameters, is also saved in _pdraws.
0005 % This routine is more efficient than metropolis_draw.m because here an
0006 % _mh file cannot be opened twice.
0007 %
0008 % INPUTS
0009 %   o SampleSize     [integer]  Size of the sample to build.
0010 %   o drsize         [double]   structure dr is drsize megaoctets.
0011 %
0012 % OUTPUTS
0013 %   o SampleAddress  [integer]  A (SampleSize*4) array, each line specifies the
0014 %                               location of a posterior draw:
0015 %                                  Column 2 --> Chain number
0016 %                                  Column 3 --> (mh) File number
0017 %                                  Column 4 --> (mh) line number
0018 %
0019 % SPECIAL REQUIREMENTS
0020 %   None.
0021 %
0022 
0023 % Copyright (C) 2006-2011 Dynare Team
0024 %
0025 % This file is part of Dynare.
0026 %
0027 % Dynare is free software: you can redistribute it and/or modify
0028 % it under the terms of the GNU General Public License as published by
0029 % the Free Software Foundation, either version 3 of the License, or
0030 % (at your option) any later version.
0031 %
0032 % Dynare is distributed in the hope that it will be useful,
0033 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0034 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0035 % GNU General Public License for more details.
0036 %
0037 % You should have received a copy of the GNU General Public License
0038 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0039 
0040 global M_ options_ estim_params_ oo_
0041 
0042 % Number of parameters:
0043 npar = estim_params_.nvx;
0044 npar = npar + estim_params_.nvn;
0045 npar = npar + estim_params_.ncx;
0046 npar = npar + estim_params_.ncn;
0047 npar = npar + estim_params_.np;
0048 
0049 % Select one task:
0050 switch nargin
0051   case 1
0052     info = 0;
0053   case 2
0054     MAX_mega_bytes = 10;% Should be an option...
0055     if drsize>0
0056         info=2;
0057     else
0058         info=1;
0059     end
0060     drawsize = drsize+npar*8/1048576;
0061   otherwise
0062     error(['selec_posterior_draws:: Unexpected number of input arguments!'])
0063 end
0064 
0065 % Get informations about the mcmc:
0066 MhDirectoryName = CheckPath('metropolis',M_.dname);
0067 fname = [ MhDirectoryName '/' M_.fname];
0068 load([ fname '_mh_history.mat']);
0069 FirstMhFile = record.KeepedDraws.FirstMhFile;
0070 FirstLine = record.KeepedDraws.FirstLine;
0071 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
0072 LastMhFile = TotalNumberOfMhFiles;
0073 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
0074 NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
0075 MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
0076 mh_nblck = options_.mh_nblck;
0077 
0078 % Randomly select draws in the posterior distribution:
0079 SampleAddress = zeros(SampleSize,4);
0080 for i = 1:SampleSize
0081     ChainNumber = ceil(rand*mh_nblck);
0082     DrawNumber  = ceil(rand*NumberOfDraws);
0083     SampleAddress(i,1) = DrawNumber;
0084     SampleAddress(i,2) = ChainNumber;
0085     if DrawNumber <= MAX_nruns-FirstLine+1
0086         MhFileNumber = FirstMhFile;
0087         MhLineNumber = FirstLine+DrawNumber-1;
0088     else
0089         DrawNumber  = DrawNumber-(MAX_nruns-FirstLine+1);
0090         MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
0091         MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns;
0092     end
0093     SampleAddress(i,3) = MhFileNumber;
0094     SampleAddress(i,4) = MhLineNumber;
0095 end
0096 SampleAddress = sortrows(SampleAddress,[3 2]);
0097 
0098 % Selected draws in the posterior distribution, and if drsize>0
0099 % reduced form solutions, are saved on disk.
0100 if info
0101     if  SampleSize*drawsize <= MAX_mega_bytes% The posterior draws are saved in one file.
0102         pdraws = cell(SampleSize,info);
0103         old_mhfile = 0;
0104         old_mhblck = 0;
0105         for i = 1:SampleSize
0106             mhfile = SampleAddress(i,3);
0107             mhblck = SampleAddress(i,2);
0108             if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
0109                 load([fname '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
0110             end
0111             pdraws(i,1) = {x2(SampleAddress(i,4),:)};
0112             if info-1
0113                 set_parameters(pdraws{i,1});
0114                 [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
0115                 pdraws(i,2) = { dr };
0116             end
0117             old_mhfile = mhfile;
0118             old_mhblck = mhblck;
0119         end
0120         clear('x2')
0121         save([fname '_posterior_draws1.mat'],'pdraws')
0122     else% The posterior draws are saved in xx files.
0123         NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize);
0124         NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes);
0125         NumberOfLines = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile;
0126         linee = 0;
0127         fnum  = 1;
0128         pdraws = cell(NumberOfDrawsPerFile,info);
0129         old_mhfile = 0;
0130         old_mhblck = 0;
0131         for i=1:SampleSize
0132             linee = linee+1;
0133             mhfile = SampleAddress(i,3);
0134             mhblck = SampleAddress(i,2);
0135             if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
0136                 load([fname '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
0137             end
0138             pdraws(linee,1) = {x2(SampleAddress(i,4),:)};
0139             if info-1
0140                 set_parameters(pdraws{linee,1});
0141                 [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
0142                 pdraws(linee,2) = { dr };
0143             end
0144             old_mhfile = mhfile;
0145             old_mhblck = mhblck;
0146             if fnum < NumberOfFiles && linee == NumberOfDrawsPerFile
0147                 linee = 0;
0148                 save([fname '_posterior_draws' num2str(fnum) '.mat'],'pdraws')
0149                 fnum = fnum+1;
0150                 if fnum < NumberOfFiles
0151                     pdraws = cell(NumberOfDrawsPerFile,info);
0152                 else
0153                     pdraws = cell(NumberOfLines,info);
0154                 end
0155             end
0156         end
0157         save([fname '_posterior_draws' num2str(fnum) '.mat'],'pdraws')
0158     end
0159 end