CutSample

 function CutSample()
 Takes a subset from metropolis

 INPUTS
   options_         [structure]
   estim_params_    [structure]
   M_               [structure]

 OUTPUTS
    none

 SPECIAL REQUIREMENTS
    none

0001 function CutSample(M_, options_, estim_params_)
0002 
0003 % function CutSample()
0004 % Takes a subset from metropolis
0005 %
0006 % INPUTS
0007 %   options_         [structure]
0008 %   estim_params_    [structure]
0009 %   M_               [structure]
0010 %
0011 % OUTPUTS
0012 %    none
0013 %
0014 % SPECIAL REQUIREMENTS
0015 %    none
0016 
0017 % Copyright (C) 2005-2011 Dynare Team
0018 %
0019 % This file is part of Dynare.
0020 %
0021 % Dynare is free software: you can redistribute it and/or modify
0022 % it under the terms of the GNU General Public License as published by
0023 % the Free Software Foundation, either version 3 of the License, or
0024 % (at your option) any later version.
0025 %
0026 % Dynare is distributed in the hope that it will be useful,
0027 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0028 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0029 % GNU General Public License for more details.
0030 %
0031 % You should have received a copy of the GNU General Public License
0032 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0033 
0034 npar = estim_params_.np+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nvx;
0035 
0036 DirectoryName = CheckPath('metropolis',M_.dname);
0037 file = dir([ DirectoryName ,filesep,  M_.fname '_mh_history.mat']);
0038 files = dir([ DirectoryName ,filesep, M_.fname '_mh*.mat' ]);
0039 if ~length(files)
0040     disp('MH:: FAILURE! there is no MH file to load here!')
0041     return
0042 end
0043 if ~length(file)
0044     disp('MH:: FAILURE! there is no MH-history file!')
0045     return
0046 else
0047     load([ DirectoryName '/'  M_.fname '_mh_history.mat'])
0048 end
0049 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
0050 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
0051 MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
0052 FirstDraw = max(1,floor(options_.mh_drop*TotalNumberOfMhDraws));
0053 FirstMhFile = ceil(FirstDraw/MAX_nruns);
0054 FirstLine = FirstDraw-(FirstMhFile-1)*MAX_nruns;
0055 record.KeepedDraws.FirstMhFile = FirstMhFile;
0056 record.KeepedDraws.FirstLine = FirstLine;
0057 if (TotalNumberOfMhFiles-1)-(FirstMhFile+1)+1 > 0
0058     record.KeepedDraws.Distribution = [ MAX_nruns-FirstLine+1 ; ...
0059                         ones((TotalNumberOfMhFiles-1)-(FirstMhFile+1)+1,1)*MAX_nruns ; ...
0060                         record.MhDraws(end,3) ];
0061 elseif TotalNumberOfMhFiles == 1
0062     record.KeepedDraws.Distribution = [];
0063 elseif TotalNumberOfMhFiles == 2 && FirstMhFile > 1
0064     record.KeepedDraws.Distribution = [MAX_nruns-FirstLine+1 ; record.MhDraws(end,3)];  
0065 end
0066 save([DirectoryName '/' M_.fname '_mh_history.mat'],'record');
0067 fprintf('MH: Total number of Mh draws: %d.\n',TotalNumberOfMhDraws);
0068 fprintf('MH: Total number of generated Mh files: %d.\n',TotalNumberOfMhFiles);
0069 fprintf('MH: I''ll use mh-files %d to %d.\n',FirstMhFile,TotalNumberOfMhFiles);
0070 fprintf('MH: In mh-file number %d i''ll start at line %d.\n',FirstMhFile,FirstLine);
0071 fprintf('MH: Finally I keep %d draws.\n',TotalNumberOfMhDraws-FirstDraw);
0072 disp(' ');