marginal_density

 function marginal = marginal_density()
 Computes the marginal density

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

 OUTPUTS
   marginal:        [double]     marginal density (modified harmonic mean)
   oo_              [structure]

 SPECIAL REQUIREMENTS
    none

0001 function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_)
0002 % function marginal = marginal_density()
0003 % Computes the marginal density
0004 %
0005 % INPUTS
0006 %   options_         [structure]
0007 %   estim_params_    [structure]
0008 %   M_               [structure]
0009 %   oo_              [structure]
0010 %
0011 % OUTPUTS
0012 %   marginal:        [double]     marginal density (modified harmonic mean)
0013 %   oo_              [structure]
0014 %
0015 % SPECIAL REQUIREMENTS
0016 %    none
0017 
0018 % Copyright (C) 2005-2011 Dynare Team
0019 %
0020 % This file is part of Dynare.
0021 %
0022 % Dynare is free software: you can redistribute it and/or modify
0023 % it under the terms of the GNU General Public License as published by
0024 % the Free Software Foundation, either version 3 of the License, or
0025 % (at your option) any later version.
0026 %
0027 % Dynare is distributed in the hope that it will be useful,
0028 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0029 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0030 % GNU General Public License for more details.
0031 %
0032 % You should have received a copy of the GNU General Public License
0033 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0034 
0035 
0036 npar = estim_params_.np+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nvx;
0037 nblck = options_.mh_nblck;
0038 
0039 MhDirectoryName = CheckPath('metropolis',M_.dname);
0040 load([ MhDirectoryName '/'  M_.fname '_mh_history.mat'])
0041 
0042 FirstMhFile = record.KeepedDraws.FirstMhFile;
0043 FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;
0044 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
0045 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
0046 MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
0047 TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws);
0048 
0049 fprintf('MH: I''m computing the posterior mean and covariance... ');
0050 [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix();
0051 
0052 MU = transpose(posterior_mean);
0053 SIGMA = posterior_covariance;
0054 lpost_mode = posterior_kernel_at_the_mode;
0055 xparam1 = posterior_mean;
0056 hh = inv(SIGMA);
0057 fprintf(' Done!\n');
0058 
0059 %% save the posterior mean and the inverse of the covariance matrix
0060 %% (usefull if the user wants to perform some computations using
0061 %% the posterior mean instead of the posterior mode ==> ).
0062 save([M_.fname '_mean.mat'],'xparam1','hh','SIGMA');
0063 
0064 disp(' ');
0065 disp('MH: I''m computing the posterior log marginale density (modified harmonic mean)... ');
0066 detSIGMA = det(SIGMA);
0067 invSIGMA = inv(SIGMA);
0068 marginal = zeros(9,2);
0069 linee = 0;
0070 check_coverage = 1;
0071 increase = 1;
0072 while check_coverage
0073     for p = 0.1:0.1:0.9;
0074         critval = chi2inv(p,npar);
0075         ifil = FirstLine;
0076         tmp = 0;
0077         for n = FirstMhFile:TotalNumberOfMhFiles
0078             for b=1:nblck
0079                 load([ MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2','logpo2');
0080                 EndOfFile = size(x2,1);
0081                 for i = ifil:EndOfFile
0082                     deviation  = (x2(i,:)-MU)*invSIGMA*(x2(i,:)-MU)';
0083                     if deviation <= critval
0084                         lftheta = -log(p)-(npar*log(2*pi)+log(detSIGMA)+deviation)/2;
0085                         tmp = tmp + exp(lftheta - logpo2(i) + lpost_mode);
0086                     end
0087                 end
0088             end
0089             ifil = 1;
0090         end
0091         linee = linee + 1;
0092         warning_old_state = warning;
0093         warning off;
0094         marginal(linee,:) = [p, lpost_mode-log(tmp/((TotalNumberOfMhDraws-TODROP)*nblck))];
0095         warning(warning_old_state);
0096     end
0097     if abs((marginal(9,2)-marginal(1,2))/marginal(9,2)) > 0.01 || isinf(marginal(1,2))
0098         if increase == 1
0099             disp('MH: The support of the weighting density function is not large enough...')
0100             disp('MH: I increase the variance of this distribution.')
0101             increase = 1.2*increase;
0102             invSIGMA = inv(SIGMA*increase);
0103             detSIGMA = det(SIGMA*increase);
0104             linee    = 0;   
0105         else
0106             disp('MH: Let me try again.')
0107             increase = 1.2*increase;
0108             invSIGMA = inv(SIGMA*increase);
0109             detSIGMA = det(SIGMA*increase);
0110             linee    = 0;
0111             if increase > 20
0112                 check_coverage = 0;
0113                 clear invSIGMA detSIGMA increase;
0114                 disp('MH: There''s probably a problem with the modified harmonic mean estimator.')    
0115             end
0116         end  
0117     else
0118         check_coverage = 0;
0119         clear invSIGMA detSIGMA increase;
0120         disp('MH: Modified harmonic mean estimator, done!')
0121     end
0122 end
0123 
0124 oo_.MarginalDensity.ModifiedHarmonicMean = mean(marginal(:,2));