From 6d08e650711e59ea741bc5206bb3f91d7d9404b8 Mon Sep 17 00:00:00 2001 From: adjemian Date: Thu, 18 Sep 2008 16:28:57 +0000 Subject: [PATCH] v4.1: Code factorization and cosmetic changes, the marginal density use a recursive approach to compute the posterior mean and covariance. git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2083 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/compute_mh_covariance_matrix.m | 48 +++++++++++++++------------ matlab/marginal_density.m | 37 +++++---------------- 2 files changed, 34 insertions(+), 51 deletions(-) diff --git a/matlab/compute_mh_covariance_matrix.m b/matlab/compute_mh_covariance_matrix.m index bb458c872..25f759b62 100644 --- a/matlab/compute_mh_covariance_matrix.m +++ b/matlab/compute_mh_covariance_matrix.m @@ -1,4 +1,4 @@ -function [m0,s0] = compute_mh_covariance_matrix() +function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix() % function [m0,s0] = compute_mh_covariance_matrix() % Estimation of the posterior covariance matrix and expectation. @@ -7,9 +7,10 @@ function [m0,s0] = compute_mh_covariance_matrix() % None. % % OUTPUTS -% o m0 [double] (n*1) vector, posterior expectation of the parameters. -% o s0 [double] (n*n) matrix, posterior covariance of the parameters -% (computed from previous metropolis hastings). +% o posterior_mean [double] (n*1) vector, posterior expectation of the parameters. +% o posterior_covariance [double] (n*n) matrix, posterior covariance of the parameters (computed from previous metropolis hastings). +% o posterior_mode [double] (n*1) vector, posterior mode of the parameters. +% o posterior_kernel_at_the_mode [double] scalar. % % SPECIAL REQUIREMENTS % None. @@ -33,6 +34,7 @@ function [m0,s0] = compute_mh_covariance_matrix() global M_ options_ estim_params_ + n = estim_params_.np + ... estim_params_.nvn+ ... estim_params_.ncx+ ... @@ -44,30 +46,32 @@ MhDirectoryName = CheckPath('metropolis'); load([ MhDirectoryName '/' M_.fname '_mh_history.mat']) FirstMhFile = record.KeepedDraws.FirstMhFile; -FirstLine = record.KeepedDraws.FirstLine; +FirstLine = record.KeepedDraws.FirstLine; TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); -params = zeros(1,n); -oldlogpo2 = -Inf; +posterior_kernel_at_the_mode = -Inf; +posterior_mean = zeros(n,1); +posterior_mode = NaN(n,1); +posterior_covariance = zeros(n,n); offset = 0; -m0 = zeros(n,1); -s0 = zeros(n,n); -for n = FirstMhFile:TotalNumberOfMhFiles - for b = 1:nblck - load([ MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2','logpo2'); - [tmp,idx] = max(logpo2); - if tmp>oldlogpo2 - oldlogpo2 = tmp; - params = x2(idx,:); +for b=1:nblck + first_line = FirstLine; + for n = FirstMhFile:TotalNumberOfMhFiles + %for b = 1:nblck + load([ MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2','logpo2'); + [tmp,idx] = max(logpo2); + if tmp>posterior_kernel_at_the_mode + posterior_kernel_at_the_mode = tmp; + posterior_mode = x2(idx,:); + end + [posterior_mean,posterior_covariance,offset] = recursive_moments(posterior_mean,posterior_covariance,x2(first_line:end,:),offset); + first_line = 1; end - [m0,s0,offset] = recursive_moments(m0,s0,x2(FirstLine,:),offset); - end - FirstLine = 1; end -xparam1 = params'; -hh = inv(s0); -fval = oldlogpo2; +xparam1 = posterior_mode'; +hh = inv(posterior_covariance); +fval = posterior_kernel_at_the_mode; save([M_.fname '_mh_mode.mat'],'xparam1','hh','fval'); \ No newline at end of file diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m index a7074b997..b57a08fdc 100644 --- a/matlab/marginal_density.m +++ b/matlab/marginal_density.m @@ -46,42 +46,21 @@ TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws); -MU = zeros(1,npar); -SIGMA = zeros(npar,npar); -lpost_mode = -Inf; +fprintf('MH: I''m computing the posterior mean and covariance... '); +[posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(); -fprintf('MH: I''m computing the posterior mean... '); -for n = FirstMhFile:TotalNumberOfMhFiles - for b = 1:nblck - load([ MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2','logpo2'); - MU = MU + sum(x2(ifil:end,:)); - lpost_mode = max(lpost_mode,max(logpo2)); - end - ifil = 1; -end -MU = MU/((TotalNumberOfMhDraws-TODROP)*nblck); -xparam1 = MU'; -MU1 = repmat(MU,MAX_nruns,1); -%% lpost_mode is the value of the log posterior kernel at the mode. -fprintf(' Done!\n'); -fprintf('MH: I''m computing the posterior covariance matrix... '); -ifil = FirstLine; -for n = FirstMhFile:TotalNumberOfMhFiles - for b = 1:nblck - load([MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2'); - x = x2(ifil:end,:)-MU1(1:size(x2(ifil:end,:),1),:); - SIGMA = SIGMA + x'*x; - end - ifil = 1; -end -SIGMA = SIGMA/((TotalNumberOfMhDraws-TODROP)*nblck);%<=== Variance of the parameters (ok!) +MU = transpose(posterior_mean); +SIGMA = posterior_covariance; +lpost_mode = posterior_kernel_at_the_mode; +xparam1 = posterior_mean; hh = inv(SIGMA); fprintf(' Done!\n'); + %% save the posterior mean and the inverse of the covariance matrix %% (usefull if the user wants to perform some computations using %% the posterior mean instead of the posterior mode ==> ). save([M_.fname '_mean.mat'],'xparam1','hh','SIGMA'); -%% end%Save. + disp(' '); disp('MH: I''m computing the posterior log marginale density (modified harmonic mean)... '); detSIGMA = det(SIGMA);