2009-12-16 18:17:34 +01:00
function myoutput = McMCDiagnostics_core ( myinputs,fpar,npar,whoiam, ThisMatlab)
2014-10-09 10:57:28 +02:00
% function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
% Computes the Brooks/Gelman (1998) convergence diagnostics, both the
% parameteric and the non-parameteric versions
2017-05-16 15:10:20 +02:00
%
2010-05-31 11:47:30 +02:00
% PARALLEL CONTEXT
% Core functionality for MCMC Diagnostics, which can be parallelized.
2016-05-19 17:13:39 +02:00
% See also the comment in posterior_sampler_core.m funtion.
2017-05-16 15:10:20 +02:00
%
%
2011-02-02 14:13:11 +01:00
% INPUTS
2016-05-19 17:13:39 +02:00
% See See the comment in posterior_sampler_core.m funtion.
2010-05-31 11:47:30 +02:00
% OUTPUTS
% o myoutput [struc]
2014-10-09 10:57:28 +02:00
% Contains:
% - UDIAG [by 6] double 1st column: length of total sequence interval
% 2nd column: sum of length of within sequence intervals; used to compute mean length of within sequence intervals
% 3nd column: within sequence variance
% 4nd column: sum of within sequence variances; used to compute mean within sequence variances
% 5nd column: within sequence kurtosis
% 6nd column: sum of within sequence kurtoses; used to compute mean within sequence kurtoses
% Averaging to compute mean moments is done in
% McMCDiagnostics
2010-05-31 11:47:30 +02:00
%
2011-02-02 14:13:11 +01:00
% ALGORITHM
2014-10-09 10:57:28 +02:00
% Computes part of the convergence diagnostics, the rest is computed in McMCDiagnostics.m .
2017-05-16 15:10:20 +02:00
% The methodology and terminology is based on: Brooks/Gelman (1998): General
% Methods for Monitoring Convergence of Iterative Simulations, Journal of Computational
2014-10-09 10:57:28 +02:00
% and Graphical Statistics, Volume 7, Number 4, Pages 434-455
2017-05-16 15:10:20 +02:00
%
2010-05-31 11:47:30 +02:00
%
% SPECIAL REQUIREMENTS.
% None.
2017-05-16 14:11:15 +02:00
% Copyright (C) 2006-2017 Dynare Team
2009-12-16 18:17:34 +01:00
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
2017-05-16 12:42:01 +02:00
if nargin < 4
2009-12-16 18:17:34 +01:00
whoiam = 0 ;
end
2010-05-31 11:47:30 +02:00
% Reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
2013-11-22 11:28:24 +01:00
MetropolisFolder = myinputs . MetropolisFolder ; %myinputs.MetropolisFolder;
2010-05-31 11:47:30 +02:00
nblck = myinputs . nblck ;
NumberOfMcFilesPerBlock = myinputs . NumberOfMcFilesPerBlock ;
Origin = myinputs . Origin ;
StepSize = myinputs . StepSize ;
mh_drop = myinputs . mh_drop ;
NumberOfDraws = myinputs . NumberOfDraws ;
NumberOfLines = myinputs . NumberOfLines ;
time = myinputs . time ;
M_ = myinputs . M_ ;
if whoiam
Parallel = myinputs . Parallel ;
end
2017-05-16 12:42:01 +02:00
if ~ exist ( ' MetropolisFolder' )
2013-11-22 11:28:24 +01:00
MetropolisFolder = CheckPath ( ' metropolis' , M_ . dname ) ;
2009-12-16 18:17:34 +01:00
end
2014-10-09 10:57:28 +02:00
ALPHA = 0.2 ; % percentile for non-parametric statistic
2009-12-16 18:17:34 +01:00
tmp = zeros ( NumberOfDraws * nblck , 3 ) ;
UDIAG = zeros ( NumberOfLines , 6 , npar - fpar + 1 ) ;
if whoiam
waitbarString = [ ' Please wait... McMCDiagnostics (' int2str ( fpar ) ' of' int2str ( npar ) ' )...' ] ;
2017-05-16 12:42:01 +02:00
if Parallel ( ThisMatlab ) . Local
2009-12-16 18:17:34 +01:00
waitbarTitle = [ ' Local ' ] ;
else
2010-10-22 11:27:26 +02:00
waitbarTitle = [ Parallel ( ThisMatlab ) . ComputerName ] ;
2011-02-02 14:13:11 +01:00
end
fMessageStatus ( 0 , whoiam , waitbarString , waitbarTitle , Parallel ( ThisMatlab ) ) ;
2009-12-16 18:17:34 +01:00
end
2017-05-16 12:42:01 +02:00
for j = fpar : npar
2013-11-04 10:54:45 +01:00
if isoctave
2017-05-16 12:42:01 +02:00
if ( whoiam == 0 )
2011-02-02 14:13:11 +01:00
printf ( ' Parameter %d... ' , j ) ;
end
else
fprintf ( ' Parameter %d... ' , j ) ;
end
2014-10-09 10:57:28 +02:00
for b = 1 : nblck %load draws from different chains into 1 matrix
2009-12-16 18:17:34 +01:00
startline = 0 ;
for n = 1 : NumberOfMcFilesPerBlock
2013-11-22 11:28:24 +01:00
load ( [ MetropolisFolder ' /' M_ . fname ' _mh' , int2str ( n ) , ' _blck' int2str ( b ) ' .mat' ] , ' x2' ) ;
2009-12-16 18:17:34 +01:00
nx2 = size ( x2 , 1 ) ;
tmp ( ( b - 1 ) * NumberOfDraws + startline + ( 1 : nx2 ) , 1 ) = x2 ( : , j ) ;
startline = startline + nx2 ;
end
end
2014-10-09 10:57:28 +02:00
tmp ( : , 2 ) = kron ( transpose ( 1 : nblck ) , ones ( NumberOfDraws , 1 ) ) ; %add info about chain associated with draw into 2nd column
tmp ( : , 3 ) = kron ( ones ( nblck , 1 ) , time ' ) ; %add timeline for draws to third column
tmp = sortrows ( tmp , 1 ) ; %sort draws according to size for non-parametric percentile computation
window_iter = 0 ;
for iter = Origin : StepSize : NumberOfDraws %begin of window
window_iter = window_iter + 1 ;
linea = ceil ( mh_drop * iter ) ; %compute first non-discarded draw; drops fraction of sample at each iteration for computational efficiency, see Brooks/Gelman (1998), p.438
n = iter - linea + 1 ; %number of draws from each block in current batch
cinf = round ( n * ALPHA / 2 ) ; %lower bound for alpha percentile of within series
csup = round ( n * ( 1 - ALPHA / 2 ) ) ; %upper bound for alpha percentile of within series
CINF = round ( nblck * n * ALPHA / 2 ) ; %lower bound for alpha percentile of pooled series
CSUP = round ( nblck * n * ( 1 - ALPHA / 2 ) ) ; %upper bound for alpha percentile of pooled series
temp = tmp ( find ( ( tmp ( : , 3 ) > = linea ) & ( tmp ( : , 3 ) < = iter ) ) , 1 : 2 ) ; %extract pooled draws in current batch
UDIAG ( window_iter , 1 , j - fpar + 1 ) = temp ( CSUP , 1 ) - temp ( CINF , 1 ) ; %length of total sequence interval
pooled_mean = mean ( temp ( : , 1 ) ) ; % Pooled mean.
UDIAG ( window_iter , 3 , j - fpar + 1 ) = sum ( ( temp ( : , 1 ) - pooled_mean ) .^ 2 ) / ( nblck * n - 1 ) ; %within sequence variance
UDIAG ( window_iter , 5 , j - fpar + 1 ) = sum ( abs ( temp ( : , 1 ) - pooled_mean ) .^ 3 ) / ( nblck * n - 1 ) ; %within sequence third moment
2009-12-16 18:17:34 +01:00
for i = 1 : nblck
pmet = temp ( find ( temp ( : , 2 ) == i ) ) ;
2014-10-09 10:57:28 +02:00
UDIAG ( window_iter , 2 , j - fpar + 1 ) = UDIAG ( window_iter , 2 , j - fpar + 1 ) + pmet ( csup , 1 ) - pmet ( cinf , 1 ) ; %sum of length of within sequence intervals; used to compute mean length of within sequence intervals
within_mean = mean ( pmet , 1 ) ; %% Within mean in current chain.
UDIAG ( window_iter , 4 , j - fpar + 1 ) = UDIAG ( window_iter , 4 , j - fpar + 1 ) + sum ( ( pmet ( : , 1 ) - within_mean ) .^ 2 ) / ( n - 1 ) ; %sum of within sequence variances; used to compute mean within sequence variances
UDIAG ( window_iter , 6 , j - fpar + 1 ) = UDIAG ( window_iter , 6 , j - fpar + 1 ) + sum ( abs ( pmet ( : , 1 ) - within_mean ) .^ 3 ) / ( n - 1 ) ; %sum of within sequence kurtoses; used to compute mean within sequence kurtoses
2009-12-16 18:17:34 +01:00
end
end
2013-11-04 10:54:45 +01:00
if isoctave
2017-05-16 12:42:01 +02:00
if ( whoiam == 0 )
2011-02-02 14:13:11 +01:00
printf ( ' Done! \n' ) ;
end
else
fprintf ( ' Done! \n' ) ;
end
2017-05-16 12:42:01 +02:00
if whoiam
2009-12-16 18:17:34 +01:00
waitbarString = [ ' Parameter ' int2str ( j ) ' /' int2str ( npar ) ' done.' ] ;
2010-10-22 11:27:26 +02:00
fMessageStatus ( ( j - fpar + 1 ) / ( npar - fpar + 1 ) , whoiam , waitbarString , waitbarTitle , Parallel ( ThisMatlab ) )
2009-12-16 18:17:34 +01:00
end
end
2009-05-27 18:30:30 +02:00
myoutput . UDIAG = UDIAG ;