diff --git a/matlab/McMCDiagnostics_core.m b/matlab/McMCDiagnostics_core.m index ad58c47e9..755d568b2 100644 --- a/matlab/McMCDiagnostics_core.m +++ b/matlab/McMCDiagnostics_core.m @@ -2,17 +2,17 @@ function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab) % PARALLEL CONTEXT % Core functionality for MCMC Diagnostics, which can be parallelized. % See also the comment in random_walk_metropolis_hastings_core.m funtion. - -% INPUTS + +% INPUTS % See See the comment in random_walk_metropolis_hastings_core.m funtion. % OUTPUTS % o myoutput [struc] % Contained UDIAG. % -% ALGORITHM -% Portion of McMCDiagnostics.m function. +% ALGORITHM +% Portion of McMCDiagnostics.m function. % % SPECIAL REQUIREMENTS. % None. @@ -59,7 +59,7 @@ if ~exist('MhDirectoryName'), MhDirectoryName = CheckPath('metropolis'); end -ALPHA = 0.2; % increase too much with the number of simulations. +ALPHA = 0.2; % increase too much with the number of simulations. tmp = zeros(NumberOfDraws*nblck,3); UDIAG = zeros(NumberOfLines,6,npar-fpar+1); @@ -69,30 +69,36 @@ if whoiam waitbarTitle=['Local ']; else waitbarTitle=[Parallel(ThisMatlab).ComputerName]; - end - fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab)); + end + fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab)); end for j=fpar:npar, - fprintf(' Parameter %d... ',j); + if exist('OCTAVE_VERSION'), + if (whoiam==0), + printf(' Parameter %d... ',j); + end + else + fprintf(' Parameter %d... ',j); + end for b = 1:nblck startline = 0; for n = 1:NumberOfMcFilesPerBlock %load([MhDirectoryName '/' mcfiles(n,1,b).name],'x2'); load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) ... - '.mat'],'x2'); + '.mat'],'x2'); nx2 = size(x2,1); tmp((b-1)*NumberOfDraws+startline+(1:nx2),1) = x2(:,j); % clear x2; startline = startline + nx2; end -% $$$ %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'x2'); -% $$$ load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'x2'); -% $$$ tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = x2(:,j); -% $$$ clear x2; -% $$$ startline = startline + LastLineNumber; + % $$$ %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'x2'); + % $$$ load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'x2'); + % $$$ tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = x2(:,j); + % $$$ clear x2; + % $$$ startline = startline + LastLineNumber; end tmp(:,2) = kron(transpose(1:nblck),ones(NumberOfDraws,1)); - tmp(:,3) = kron(ones(nblck,1),time'); + tmp(:,3) = kron(ones(nblck,1),time'); tmp = sortrows(tmp,1); ligne = 0; for iter = Origin:StepSize:NumberOfDraws @@ -111,13 +117,19 @@ for j=fpar:npar, for i=1:nblck pmet = temp(find(temp(:,2)==i)); UDIAG(ligne,2,j-fpar+1) = UDIAG(ligne,2,j-fpar+1) + pmet(csup,1)-pmet(cinf,1); - moyenne = mean(pmet,1); %% Within mean. + moyenne = mean(pmet,1); %% Within mean. UDIAG(ligne,4,j-fpar+1) = UDIAG(ligne,4,j-fpar+1) + sum((pmet(:,1)-moyenne).^2)/(n-1); UDIAG(ligne,6,j-fpar+1) = UDIAG(ligne,6,j-fpar+1) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1); end end - fprintf('Done! \n'); - if whoiam, + if exist('OCTAVE_VERSION'), + if (whoiam==0), + printf('Done! \n'); + end + else + fprintf('Done! \n'); + end + if whoiam, waitbarString = [ 'Parameter ' int2str(j) '/' int2str(npar) ' done.']; fMessageStatus((j-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab)) end diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m index 5d4732563..cb8fa57c3 100644 --- a/matlab/PosteriorIRF_core1.m +++ b/matlab/PosteriorIRF_core1.m @@ -5,10 +5,10 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,npar,whoiam, ThisMatlab) % that running in parallel a 'for' cycle, this function run in parallel a % 'while' loop! The parallelization of 'while' loop (when possible) is a more % sophisticated procedure. -% +% % See also the comment in random_walk_metropolis_hastings_core.m funtion. % -% INPUTS +% INPUTS % See the comment in random_walk_metropolis_hastings_core.m funtion. % % OUTPUTS @@ -16,8 +16,8 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,npar,whoiam, ThisMatlab) % Contained: % OutputFileName_dsge, OutputFileName_param and OutputFileName_bvardsge. % -% ALGORITHM -% Portion of PosteriorIRF.m function. Specifically the 'while' cycle. +% ALGORITHM +% Portion of PosteriorIRF.m function. Specifically the 'while' cycle. % % SPECIAL REQUIREMENTS. % None. @@ -109,10 +109,10 @@ else h = waitbar(0,'Bayesian (prior) IRFs...'); end end - + end -OutputFileName_bvardsge = {}; +OutputFileName_bvardsge = {}; OutputFileName_dsge = {}; OutputFileName_param = {}; @@ -121,25 +121,25 @@ fpar0=fpar; fpar = fpar-1; if whoiam - ifil2=ifil2(whoiam); - NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge(whoiam); - NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar(whoiam); + ifil2=ifil2(whoiam); + NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge(whoiam); + NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar(whoiam); end % Parallel 'while' very good!!! -while fpar 1e-12 + if max(y(IndxVariables(j),:)) - min(y(IndxVariables(j),:)) > 1e-12 stock_irf_dsge(:,j,i,irun) = transpose(y(IndxVariables(j),:)); end end @@ -179,21 +179,21 @@ while fpar1.000000001); end - % Get the mean - mu = zeros(1,nvobs); + % Get the mean + mu = zeros(1,nvobs); % Get rotation if dsge_prior_weight > 0 Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e); @@ -213,7 +213,7 @@ while fpar 1e-10 + if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10 tmp_dsgevar(:,j) = (irfs(:,j)); end end @@ -222,12 +222,12 @@ while fpar1) & OpenOldFile(b) load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ... - '_blck' int2str(b) '.mat']) + '_blck' int2str(b) '.mat']) x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)]; logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)]; OpenOldFile(b) = 0; @@ -124,22 +124,22 @@ for b = fblck:nblck, waitbarTitle=['Local ']; else waitbarTitle=[options_.parallel(ThisMatlab).ComputerName]; - end + end fMessageStatus(0,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab)); else, hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']); set(hh,'Name','Metropolis-Hastings'); - + end isux = 0; jsux = 0; irun = fline(b); j = 1; while j <= nruns(b) - par = feval(ProposalFun, xparam1, proposal_covariance, n); + par = feval(ProposalFun, xparam1, proposal_covariance, n); if all( par(:) > mh_bounds(:,1) ) & all( par(:) < mh_bounds(:,2) ) try - logpost = - feval(TargetFun, par(:),varargin{:}); + logpost = - feval(TargetFun, par(:),varargin{:}); catch, logpost = -inf; end @@ -152,11 +152,11 @@ for b = fblck:nblck, if (logpost > -inf) && (log(rand) < r) x2(irun,:) = par; ix2(b,:) = par; - logpo2(irun) = logpost; + logpo2(irun) = logpost; ilogpo2(b) = logpost; isux = isux + 1; jsux = jsux + 1; - else + else x2(irun,:) = ix2(b,:); logpo2(irun) = ilogpo2(b); end @@ -164,12 +164,14 @@ for b = fblck:nblck, if exist('OCTAVE_VERSION') || options_.console_mode if mod(j, 10) == 0 if exist('OCTAVE_VERSION') - printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j); + if (whoiam==0), + printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j); + end else fprintf(' MH: Computing Metropolis-Hastings (chain %d/%d): %3.f \b%% done, acception rate: %3.f \b%%\r', b, nblck, 100 * prtfrc, 100 * isux / j); end end - if mod(j,50)==0 & whoiam, + if mod(j,50)==0 & whoiam, % keyboard; waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%%%', 100 * isux/j)]; fMessageStatus(prtfrc,whoiam,waitbarString, '', options_.parallel(ThisMatlab)) @@ -177,13 +179,13 @@ for b = fblck:nblck, else if mod(j, 3)==0 & ~whoiam waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]); - elseif mod(j,50)==0 & whoiam, + elseif mod(j,50)==0 & whoiam, % keyboard; waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]; fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab)) end end - + if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2'); fidlog = fopen([MhDirectoryName '/metropolis.log'],'a'); diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m index b6a324a65..c8dbb278f 100644 --- a/matlab/prior_posterior_statistics_core.m +++ b/matlab/prior_posterior_statistics_core.m @@ -3,22 +3,22 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa % Core functionality for prior_posterior.m function, which can be parallelized. % See also the comment in random_walk_metropolis_hastings_core.m funtion. % -% INPUTS +% INPUTS % See See the comment in random_walk_metropolis_hastings_core.m funtion. % OUTPUTS % o myoutput [struc] -% Contained OutputFileName_smooth; -% _update; -% _inno; -% _error; -% _filter_step_ahead; +% Contained OutputFileName_smooth; +% _update; +% _inno; +% _error; +% _filter_step_ahead; % _param; % _forc_mean; % _forc_point % -% ALGORITHM -% Portion of prior_posterior.m function. +% ALGORITHM +% Portion of prior_posterior.m function. % This file is part of Dynare. % % SPECIAL REQUIREMENTS. @@ -111,37 +111,37 @@ else else h = waitbar(0,'Taking subdraws...'); end - + end if RemoteFlag==1, -OutputFileName_smooth = {}; -OutputFileName_update = {}; -OutputFileName_inno = {}; -OutputFileName_error = {}; -OutputFileName_filter_step_ahead = {}; -OutputFileName_param = {}; -OutputFileName_forc_mean = {}; -OutputFileName_forc_point = {}; -% OutputFileName_moments = {}; + OutputFileName_smooth = {}; + OutputFileName_update = {}; + OutputFileName_inno = {}; + OutputFileName_error = {}; + OutputFileName_filter_step_ahead = {}; + OutputFileName_param = {}; + OutputFileName_forc_mean = {}; + OutputFileName_forc_point = {}; + % OutputFileName_moments = {}; end for b=fpar:B - + % [deep, logpo] = GetOneDraw(type); % set_all_parameters(deep); % dr = resol(oo_.steady_state,0); if strcmpi(type,'prior') - + [deep, logpo] = GetOneDraw(type); - + else deep = x(b,:); logpo = logpost(b); end set_all_parameters(deep); [dr,info] = resol(oo_.steady_state,0); - + if run_smoother [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK] = ... DsgeSmoother(deep,gend,Y,data_index,missing_value); @@ -164,7 +164,7 @@ for b=fpar:B if naK stock_filter_step_ahead(:,dr.order_var,:,irun(4)) = aK(options_.filter_step_ahead,1:endo_nbr,:); end - + if horizon yyyy = alphahat(iendo,i_last_obs); yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr)); @@ -190,19 +190,19 @@ for b=fpar:B else yf1 = yf1 + repmat(SteadyState',[horizon+maxlag,1,1]); end - + stock_forcst_mean(:,:,irun(6)) = yf'; stock_forcst_point(:,:,irun(7)) = yf1'; end - + end stock_param(irun(5),:) = deep; stock_logpo(irun(5),1) = logpo; stock_ys(irun(5),:) = SteadyState'; - + irun = irun + ones(7,1); - - + + if irun(1) > MAX_nsmoo || b == B stock = stock_smooth(:,:,1:irun(1)-1); ifil(1) = ifil(1) + 1; @@ -211,72 +211,72 @@ for b=fpar:B stock = stock_update(:,:,1:irun(1)-1); save([DirectoryName '/' M_.fname '_update' int2str(ifil(1)) '.mat'],'stock'); if RemoteFlag==1, - OutputFileName_smooth = [OutputFileName_smooth; {[DirectoryName filesep], [M_.fname '_smooth' int2str(ifil(1)) '.mat']}]; - OutputFileName_update = [OutputFileName_update; {[DirectoryName filesep], [M_.fname '_update' int2str(ifil(1)) '.mat']}]; + OutputFileName_smooth = [OutputFileName_smooth; {[DirectoryName filesep], [M_.fname '_smooth' int2str(ifil(1)) '.mat']}]; + OutputFileName_update = [OutputFileName_update; {[DirectoryName filesep], [M_.fname '_update' int2str(ifil(1)) '.mat']}]; end irun(1) = 1; end - + if irun(2) > MAX_ninno || b == B stock = stock_innov(:,:,1:irun(2)-1); ifil(2) = ifil(2) + 1; save([DirectoryName '/' M_.fname '_inno' int2str(ifil(2)) '.mat'],'stock'); if RemoteFlag==1, - OutputFileName_inno = [OutputFileName_inno; {[DirectoryName filesep], [M_.fname '_inno' int2str(ifil(2)) '.mat']}]; + OutputFileName_inno = [OutputFileName_inno; {[DirectoryName filesep], [M_.fname '_inno' int2str(ifil(2)) '.mat']}]; end irun(2) = 1; end - + if nvn && (irun(3) > MAX_nerro || b == B) stock = stock_error(:,:,1:irun(3)-1); ifil(3) = ifil(3) + 1; save([DirectoryName '/' M_.fname '_error' int2str(ifil(3)) '.mat'],'stock'); if RemoteFlag==1, - OutputFileName_error = [OutputFileName_error; {[DirectoryName filesep], [M_.fname '_error' int2str(ifil(3)) '.mat']}]; + OutputFileName_error = [OutputFileName_error; {[DirectoryName filesep], [M_.fname '_error' int2str(ifil(3)) '.mat']}]; end irun(3) = 1; end - + if naK && (irun(4) > MAX_naK || b == B) stock = stock_filter_step_ahead(:,:,:,1:irun(4)-1); ifil(4) = ifil(4) + 1; save([DirectoryName '/' M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat'],'stock'); if RemoteFlag==1, - OutputFileName_filter_step_ahead = [OutputFileName_filter_step_ahead; {[DirectoryName filesep], [M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat']}]; + OutputFileName_filter_step_ahead = [OutputFileName_filter_step_ahead; {[DirectoryName filesep], [M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat']}]; end irun(4) = 1; end - + if irun(5) > MAX_nruns || b == B stock = stock_param(1:irun(5)-1,:); ifil(5) = ifil(5) + 1; save([DirectoryName '/' M_.fname '_param' int2str(ifil(5)) '.mat'],'stock','stock_logpo','stock_ys'); if RemoteFlag==1, - OutputFileName_param = [OutputFileName_param; {[DirectoryName filesep], [M_.fname '_param' int2str(ifil(5)) '.mat']}]; + OutputFileName_param = [OutputFileName_param; {[DirectoryName filesep], [M_.fname '_param' int2str(ifil(5)) '.mat']}]; end irun(5) = 1; end - + if horizon && (irun(6) > MAX_nforc1 || b == B) stock = stock_forcst_mean(:,:,1:irun(6)-1); ifil(6) = ifil(6) + 1; save([DirectoryName '/' M_.fname '_forc_mean' int2str(ifil(6)) '.mat'],'stock'); if RemoteFlag==1, - OutputFileName_forc_mean = [OutputFileName_forc_mean; {[DirectoryName filesep], [M_.fname '_forc_mean' int2str(ifil(6)) '.mat']}]; + OutputFileName_forc_mean = [OutputFileName_forc_mean; {[DirectoryName filesep], [M_.fname '_forc_mean' int2str(ifil(6)) '.mat']}]; end irun(6) = 1; end - + if horizon && (irun(7) > MAX_nforc2 || b == B) stock = stock_forcst_point(:,:,1:irun(7)-1); ifil(7) = ifil(7) + 1; save([DirectoryName '/' M_.fname '_forc_point' int2str(ifil(7)) '.mat'],'stock'); if RemoteFlag==1, - OutputFileName_forc_point = [OutputFileName_forc_point; {[DirectoryName filesep], [M_.fname '_forc_point' int2str(ifil(7)) '.mat']}]; + OutputFileName_forc_point = [OutputFileName_forc_point; {[DirectoryName filesep], [M_.fname '_forc_point' int2str(ifil(7)) '.mat']}]; end irun(7) = 1; end - + % if moments_varendo && (irun(8) > MAX_momentsno || b == B) % stock = stock_moments(1:irun(8)-1); % ifil(8) = ifil(8) + 1; @@ -286,18 +286,22 @@ for b=fpar:B % end % irun(8) = 1; % end - -% DirectoryName=TempPath; - - - if exist('OCTAVE_VERSION') - printf('Taking subdraws: %3.f%% done\r', b/B*100); + + % DirectoryName=TempPath; + + + if exist('OCTAVE_VERSION'), + if (whoiam==0), + printf('Taking subdraws: %3.f%% done\r', b/B*100); + end elseif ~whoiam, waitbar(b/B,h); end - + if whoiam, - fprintf('Done! \n'); + if ~exist('OCTAVE_VERSION') + fprintf('Done! \n'); + end waitbarString = [ 'Subdraw ' int2str(b) '/' int2str(B) ' done.']; fMessageStatus((b-fpar+1)/(B-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab)); end @@ -305,15 +309,15 @@ end myoutput.ifil=ifil; if RemoteFlag==1, -myoutput.OutputFileName = [OutputFileName_smooth; -OutputFileName_update; -OutputFileName_inno; -OutputFileName_error; -OutputFileName_filter_step_ahead; -OutputFileName_param; -OutputFileName_forc_mean; -OutputFileName_forc_point]; -% OutputFileName_moments]; + myoutput.OutputFileName = [OutputFileName_smooth; + OutputFileName_update; + OutputFileName_inno; + OutputFileName_error; + OutputFileName_filter_step_ahead; + OutputFileName_param; + OutputFileName_forc_mean; + OutputFileName_forc_point]; + % OutputFileName_moments]; end if exist('OCTAVE_VERSION') @@ -323,5 +327,5 @@ else if exist('h') close(h) end - + end diff --git a/matlab/random_walk_metropolis_hastings_core.m b/matlab/random_walk_metropolis_hastings_core.m index f903ce1be..0e3a1750c 100644 --- a/matlab/random_walk_metropolis_hastings_core.m +++ b/matlab/random_walk_metropolis_hastings_core.m @@ -3,8 +3,8 @@ function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,wh % This function contain the most computationally intensive portion of code in % random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in 'for' % cycle and are completely independent than suitable to be executed in parallel way. -% -% INPUTS +% +% INPUTS % o myimput [struc] The mandatory variables for local/remote % parallel computing obtained from random_walk_metropolis_hastings.m % function. @@ -26,8 +26,8 @@ function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,wh % NewFile; % OutputFileName % -% ALGORITHM -% Portion of Metropolis-Hastings. +% ALGORITHM +% Portion of Metropolis-Hastings. % % SPECIAL REQUIREMENTS. % None. @@ -88,16 +88,16 @@ nruns=myinputs.nruns; NewFile=myinputs.NewFile; MAX_nruns=myinputs.MAX_nruns; d=myinputs.d; -InitSizeArray=myinputs.InitSizeArray; +InitSizeArray=myinputs.InitSizeArray; record=myinputs.record; varargin=myinputs.varargin; % Necessary only for remote computing! if whoiam - Parallel=myinputs.Parallel; - % initialize persistent variables in priordens() - priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ... - bayestopt_.p3,bayestopt_.p4,1); + Parallel=myinputs.Parallel; + % initialize persistent variables in priordens() + priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ... + bayestopt_.p3,bayestopt_.p4,1); end % (re)Set the penalty @@ -138,7 +138,7 @@ for b = fblck:nblck, rand('state',record.Seeds(b).Unifor); if (options_.load_mh_file~=0) & (fline(b)>1) & OpenOldFile(b) load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ... - '_blck' int2str(b) '.mat']) + '_blck' int2str(b) '.mat']) x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)]; logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)]; OpenOldFile(b) = 0; @@ -160,14 +160,14 @@ for b = fblck:nblck, hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']); set(hh,'Name','Metropolis-Hastings'); end - + end if whoiam if options_.parallel(ThisMatlab).Local, waitbarTitle=['Local ']; else waitbarTitle=[options_.parallel(ThisMatlab).ComputerName]; - end + end prc0=(b-fblck)/(nblck-fblck+1)*(exist('OCTAVE_VERSION') || options_.console_mode); fMessageStatus(prc0,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab)); end @@ -179,7 +179,7 @@ for b = fblck:nblck, par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n); if all( par(:) > mh_bounds(:,1) ) & all( par(:) < mh_bounds(:,2) ) try - logpost = - feval(TargetFun, par(:),varargin{:}); + logpost = - feval(TargetFun, par(:),varargin{:}); catch logpost = -inf; end @@ -189,11 +189,11 @@ for b = fblck:nblck, if (logpost > -inf) && (log(rand) < logpost-ilogpo2(b)) x2(irun,:) = par; ix2(b,:) = par; - logpo2(irun) = logpost; + logpo2(irun) = logpost; ilogpo2(b) = logpost; isux = isux + 1; jsux = jsux + 1; - else + else x2(irun,:) = ix2(b,:); logpo2(irun) = ilogpo2(b); end @@ -201,14 +201,16 @@ for b = fblck:nblck, if exist('OCTAVE_VERSION') || options_.console_mode if mod(j, 10) == 0 if exist('OCTAVE_VERSION') - printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j); + if (whoiam==0) + printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j); + end else s0=repmat('\b',1,length(newString)); newString=sprintf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acceptance rate: %3.f%%', b, nblck, 100 * prtfrc, 100 * isux / j); fprintf([s0,'%s'],newString); end end - if mod(j,50)==0 & whoiam + if mod(j,50)==0 & whoiam % keyboard; waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%', 100 * isux/j)]; fMessageStatus((b-fblck)/(nblck-fblck+1)+prtfrc/(nblck-fblck+1),whoiam,waitbarString, '', options_.parallel(ThisMatlab)); @@ -216,13 +218,13 @@ for b = fblck:nblck, else if mod(j, 3)==0 & ~whoiam waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]); - elseif mod(j,50)==0 & whoiam, + elseif mod(j,50)==0 & whoiam, % keyboard; waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]; fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab)); end end - + if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2'); fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');