McMCDiagnostics_core

 PARALLEL CONTEXT
 Core functionality for MCMC Diagnostics, which can be parallelized.
 See also the comment in random_walk_metropolis_hastings_core.m funtion.

0001 function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
0002 % PARALLEL CONTEXT
0003 % Core functionality for MCMC Diagnostics, which can be parallelized.
0004 % See also the comment in random_walk_metropolis_hastings_core.m funtion.
0005 
0006 
0007 % INPUTS
0008 %   See See the comment in random_walk_metropolis_hastings_core.m funtion.
0009 
0010 % OUTPUTS
0011 % o myoutput  [struc]
0012 %  Contained UDIAG.
0013 %
0014 % ALGORITHM
0015 %   Portion of McMCDiagnostics.m function.
0016 %
0017 % SPECIAL REQUIREMENTS.
0018 %   None.
0019 
0020 % Copyright (C) 2006-2011 Dynare Team
0021 %
0022 % This file is part of Dynare.
0023 %
0024 % Dynare is free software: you can redistribute it and/or modify
0025 % it under the terms of the GNU General Public License as published by
0026 % the Free Software Foundation, either version 3 of the License, or
0027 % (at your option) any later version.
0028 %
0029 % Dynare is distributed in the hope that it will be useful,
0030 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0031 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0032 % GNU General Public License for more details.
0033 %
0034 % You should have received a copy of the GNU General Public License
0035 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0036 
0037 if nargin<4,
0038     whoiam=0;
0039 end
0040 
0041 % Reshape 'myinputs' for local computation.
0042 % In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
0043 
0044 MhDirectoryName=myinputs.MhDirectoryName;
0045 nblck=myinputs.nblck;
0046 NumberOfMcFilesPerBlock=myinputs.NumberOfMcFilesPerBlock;
0047 Origin=myinputs.Origin;
0048 StepSize=myinputs.StepSize;
0049 mh_drop=myinputs.mh_drop;
0050 NumberOfDraws=myinputs.NumberOfDraws;
0051 NumberOfLines=myinputs.NumberOfLines;
0052 time=myinputs.time;
0053 M_=myinputs.M_;
0054 
0055 if whoiam
0056     Parallel=myinputs.Parallel;
0057 end
0058 if ~exist('MhDirectoryName'),
0059     MhDirectoryName = CheckPath('metropolis',M_.dname);
0060 end
0061 
0062 ALPHA = 0.2;                                % increase too much with the number of simulations.
0063 tmp = zeros(NumberOfDraws*nblck,3);
0064 UDIAG = zeros(NumberOfLines,6,npar-fpar+1);
0065 
0066 if whoiam
0067     waitbarString = ['Please wait... McMCDiagnostics (' int2str(fpar) 'of' int2str(npar) ')...'];
0068     if Parallel(ThisMatlab).Local,
0069         waitbarTitle=['Local '];
0070     else
0071         waitbarTitle=[Parallel(ThisMatlab).ComputerName];
0072     end
0073     fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
0074 end
0075 for j=fpar:npar,
0076     if exist('OCTAVE_VERSION'),
0077         if (whoiam==0),
0078             printf('    Parameter %d...  ',j);
0079         end
0080     else
0081         fprintf('    Parameter %d...  ',j);
0082     end
0083     for b = 1:nblck
0084         startline = 0;
0085         for n = 1:NumberOfMcFilesPerBlock
0086             %load([MhDirectoryName '/' mcfiles(n,1,b).name],'x2');
0087             load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) ...
0088                   '.mat'],'x2');
0089             nx2 = size(x2,1);
0090             tmp((b-1)*NumberOfDraws+startline+(1:nx2),1) = x2(:,j);
0091             %      clear x2;
0092             startline = startline + nx2;
0093         end
0094 % $$$     %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'x2');
0095 % $$$     load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'x2');
0096 % $$$     tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = x2(:,j);
0097 % $$$     clear x2;
0098 % $$$     startline = startline + LastLineNumber;
0099     end
0100     tmp(:,2) = kron(transpose(1:nblck),ones(NumberOfDraws,1));
0101     tmp(:,3) = kron(ones(nblck,1),time');
0102     tmp = sortrows(tmp,1);
0103     ligne   = 0;
0104     for iter  = Origin:StepSize:NumberOfDraws
0105         ligne = ligne+1;
0106         linea = ceil(mh_drop*iter);
0107         n     = iter-linea+1;
0108         cinf  = round(n*ALPHA/2);
0109         csup  = round(n*(1-ALPHA/2));
0110         CINF  = round(nblck*n*ALPHA/2);
0111         CSUP  = round(nblck*n*(1-ALPHA/2));
0112         temp  = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2);
0113         UDIAG(ligne,1,j-fpar+1) = temp(CSUP,1)-temp(CINF,1);
0114         moyenne = mean(temp(:,1));%% Pooled mean.
0115         UDIAG(ligne,3,j-fpar+1) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1);
0116         UDIAG(ligne,5,j-fpar+1) = sum(abs(temp(:,1)-moyenne).^3)/(nblck*n-1);
0117         for i=1:nblck
0118             pmet = temp(find(temp(:,2)==i));
0119             UDIAG(ligne,2,j-fpar+1) = UDIAG(ligne,2,j-fpar+1) + pmet(csup,1)-pmet(cinf,1);
0120             moyenne = mean(pmet,1); %% Within mean.
0121             UDIAG(ligne,4,j-fpar+1) = UDIAG(ligne,4,j-fpar+1) + sum((pmet(:,1)-moyenne).^2)/(n-1);
0122             UDIAG(ligne,6,j-fpar+1) = UDIAG(ligne,6,j-fpar+1) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1);
0123         end
0124     end
0125     if exist('OCTAVE_VERSION'),
0126         if (whoiam==0),
0127             printf('Done! \n');
0128         end
0129     else
0130         fprintf('Done! \n');
0131     end
0132     if whoiam,
0133         waitbarString = [ 'Parameter ' int2str(j) '/' int2str(npar) ' done.'];
0134         fMessageStatus((j-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab))
0135     end
0136 end
0137 
0138 myoutput.UDIAG = UDIAG;