diff --git a/matlab/estimation/GetAllPosteriorDraws.m b/matlab/estimation/GetAllPosteriorDraws.m
index fdea14c93f56029116a2a2887379f0d34a63a3e2..f7983515541566aa27d188392889eae3b0c3e4b9 100644
--- a/matlab/estimation/GetAllPosteriorDraws.m
+++ b/matlab/estimation/GetAllPosteriorDraws.m
@@ -51,6 +51,18 @@ if issmc(options_)
else
draws = posterior.tlogpostkernel;
end
+ elseif isdsmh(options_)
+ % Load draws from the posterior distribution
+ posterior = load(sprintf('%s/dsmh/parameters_particles_final.mat', dname));
+ if column>0 || strcmp(column,'all')
+ if strcmp(column,'all')
+ draws = transpose(posterior.particles);
+ else
+ draws = transpose(posterior.particles(column,:));
+ end
+ else
+ draws = posterior.tlogpostkernel;
+ end
elseif isonline(options_)
% Load draws from the posterior distribution
posterior = load(sprintf('%s/online/parameters_particles_final.mat', dname));
diff --git a/matlab/estimation/GetPosteriorMeanVariance.m b/matlab/estimation/GetPosteriorMeanVariance.m
index 7f73899cc863c5b8d51f88a388f5fa755cc816aa..9d8638363e4381e4f06145f7fd0e6f18c4d7e19a 100644
--- a/matlab/estimation/GetPosteriorMeanVariance.m
+++ b/matlab/estimation/GetPosteriorMeanVariance.m
@@ -36,6 +36,13 @@ if issmc(options_)
mean = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
% Compute the posterior covariance
variance = (posterior.particles-mean)*(posterior.particles-mean)'/length(posterior.tlogpostkernel);
+ elseif isdsmh(options_)
+ % Load draws from the posterior distribution
+ posterior = load(sprintf('%s/dsmh/parameters_particles_final.mat', M_.dname));
+ % Compute the posterior mean
+ mean = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
+ % Compute the posterior covariance
+ variance = (posterior.particles-mean)*(posterior.particles-mean)'/length(posterior.tlogpostkernel);
elseif isonline(options_)
posterior = load(sprintf('%s/online/parameters_particles_final.mat', M_.dname));
% Compute the posterior mean
diff --git a/matlab/estimation/GetPosteriorParametersStatistics.m b/matlab/estimation/GetPosteriorParametersStatistics.m
index ebaf93666523353a3a2df5c5d54f5f594743395d..85723ea0799db1c05d84bb41b1f5309d096969f6 100644
--- a/matlab/estimation/GetPosteriorParametersStatistics.m
+++ b/matlab/estimation/GetPosteriorParametersStatistics.m
@@ -58,6 +58,9 @@ skipline()
if ishssmc(options_)
dprintf('Log data density is %f.', oo_.MarginalDensity.hssmc);
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
+elseif isdsmh(options_)
+ dprintf('Log data density is %f.', oo_.MarginalDensity.dsmh);
+ hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
elseif isonline(options_)
hpd_draws = num_draws;
elseif isdime(options_)
@@ -89,7 +92,7 @@ if np
disp(tit2)
ip = nvx+nvn+ncx+ncn+1;
for i=1:np
- if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_) || isonline(options_)
+ if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_) || isdsmh(options_) || isonline(options_)
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
name = bayestopt_.name{ip};
diff --git a/matlab/estimation/check_posterior_analysis_data.m b/matlab/estimation/check_posterior_analysis_data.m
index 6b9008d4422899c063c90cb77b20a86ab5d467cf..13d71b892541cc373ac1c3b8030cb477b8814fab 100644
--- a/matlab/estimation/check_posterior_analysis_data.m
+++ b/matlab/estimation/check_posterior_analysis_data.m
@@ -45,6 +45,8 @@ if ~issmc(options_)
else
if ishssmc(options_)
[MetropolisFolder, info] = CheckPath('hssmc',M_.dname);
+ elseif isdsmh(options_)
+ [MetropolisFolder, info] = CheckPath('dsmh',M_.dname);
elseif isonline(options_)
[MetropolisFolder, info] = CheckPath('online',M_.dname);
elseif isdime(options_)
@@ -78,6 +80,10 @@ else
% Load draws from the posterior distribution
pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
mhdate = get_date_of_a_file(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+ elseif isdsmh(options_)
+ % Load draws from the posterior distribution
+ pfiles = dir(sprintf('%s/dsmh/particles-*.mat', M_.dname));
+ mhdate = get_date_of_a_file(sprintf('%s/dsmh/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
elseif isonline(options_)
% Load draws from the posterior distribution
pfiles = dir(sprintf('%s/online/parameters_particles_final.mat', M_.dname));
diff --git a/matlab/estimation/check_posterior_sampler_options.m b/matlab/estimation/check_posterior_sampler_options.m
index fe93f77a69ec5436c5ee84ed80b28f28c25e25e3..ed9e08da286bbce38a670f1551641e8966a8794c 100644
--- a/matlab/estimation/check_posterior_sampler_options.m
+++ b/matlab/estimation/check_posterior_sampler_options.m
@@ -482,34 +482,34 @@ if init
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
- case 'proposal_distribution'
- if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
- strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
- error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
- 'rand_multivariate_student or rand_multivariate_normal as options']);
- else
- posterior_sampler_options.proposal_distribution=options_list{i,2};
- end
- case 'student_degrees_of_freedom'
- if options_list{i,2} <= 0
- error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
- else
- posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
- end
- case 'save_tmp_file'
- posterior_sampler_options.save_tmp_file = options_list{i,2};
- case 'number_of_particles'
+ case 'H'
+ posterior_sampler_options.H = options_list{i,2};
+ case 'N'
+ posterior_sampler_options.N = options_list{i,2};
+ case 'G'
+ posterior_sampler_options.G = options_list{i,2};
+ case 'K'
+ posterior_sampler_options.K = options_list{i,2};
+ case 'lambda1'
+ posterior_sampler_options.lambda1 = options_list{i,2};
+ case 'tau'
+ posterior_sampler_options.tau = options_list{i,2};
+ case 'particles'
posterior_sampler_options.particles = options_list{i,2};
+ case 'alpha0'
+ posterior_sampler_options.alpha0 = options_list{i,2};
+ case 'alpha1'
+ posterior_sampler_options.alpha1 = options_list{i,2};
otherwise
- warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
- end
+ warning(['dsmh: Unknown option (' options_list{i,1} ')!'])
end
+ end
end
-
+
options_.mode_compute = 0;
options_.cova_compute = 0;
options_.mh_replic = 0;
- options_.mh_posterior_mode_estimation = true;
+ options_.mh_posterior_mode_estimation = false;
otherwise
error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index e82e4389ddcb1e9669c851d9068a8734f33bb00d..168b72167f7b23f7741fad03ad2a891fc0bb380e 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -405,20 +405,18 @@ end
% Run SMC sampler.
%
-if isonline(options_)
- [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
- options_.posterior_sampler_options.current_options = posterior_sampler_options;
- online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_);
-elseif ishssmc(options_)
- [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
- options_.posterior_sampler_options.current_options = posterior_sampler_options;
- oo_.MarginalDensity.hssmc = hssmc(objective_function, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
-elseif isdime(options_)
+if issmc(options_)
[posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
options_.posterior_sampler_options.current_options = posterior_sampler_options;
- dime(objective_function, xparam1, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
-elseif isdsmh(options_)
- dsmh(objective_function, xparam1, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
+ if isonline(options_)
+ online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_);
+ elseif ishssmc(options_)
+ oo_.MarginalDensity.hssmc = hssmc(objective_function, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
+ elseif isdime(options_)
+ dime(objective_function, xparam1, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+ elseif isdsmh(options_)
+ oo_.MarginalDensity.dsmh = dsmh(objective_function, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
+ end
end
%
@@ -482,7 +480,7 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
oo_.lprob = trace_plot_dime(options_, M_);
end
% Estimation of the marginal density from the Mh draws:
- if ishssmc(options_) || isonline(options_) || isdime(options_) || options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
+ if isdsmh(options_) || ishssmc(options_) || isonline(options_) || isdime(options_) || options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
if ~issmc(options_)
[~, oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
end
diff --git a/matlab/estimation/get_posterior_folder_name.m b/matlab/estimation/get_posterior_folder_name.m
index fc0fd0b8e463d03b60347dc27fc6f5b1a95112d6..0036f491a24ead9b216051beb6459bacc3c2da3f 100644
--- a/matlab/estimation/get_posterior_folder_name.m
+++ b/matlab/estimation/get_posterior_folder_name.m
@@ -34,6 +34,8 @@ if ~issmc(options_)
else
if ishssmc(options_)
folder_name='hssmc';
+ elseif isdsmh(options_)
+ folder_name='dsmh';
elseif isdime(options_)
folder_name='dime';
elseif isonline(options_)
diff --git a/matlab/estimation/set_number_of_subdraws.m b/matlab/estimation/set_number_of_subdraws.m
index 0241b32050e50376e00b3f4570f7c1cbb9ec13aa..2e7ef97cc61f8a0796e00895ae8c337715614530 100644
--- a/matlab/estimation/set_number_of_subdraws.m
+++ b/matlab/estimation/set_number_of_subdraws.m
@@ -57,6 +57,11 @@ else
posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
NumberOfDraws = size(posterior.particles,2);
NumberOfDrawsPerChain=NumberOfDraws;
+ elseif isdsmh(options_)
+ % Load draws from the posterior distribution
+ posterior = load(sprintf('%s/dsmh/parameters_particles_final.mat', M_.dname));
+ NumberOfDraws = size(posterior.particles,2);
+ NumberOfDrawsPerChain=NumberOfDraws;
elseif isonline(options_)
% Load draws from the posterior distribution
posterior = load(sprintf('%s/online/parameters_particles_final.mat', M_.dname));
diff --git a/matlab/estimation/smc/dsmh.m b/matlab/estimation/smc/dsmh.m
index eb693ddb9e9a91998947186b9e8f88af67b85e2e..1059753e56ae965c2c40bb772feca9daa25854af 100644
--- a/matlab/estimation/smc/dsmh.m
+++ b/matlab/estimation/smc/dsmh.m
@@ -1,11 +1,12 @@
-function dsmh(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
-
+function mdd = dsmh(TargetFun, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
+% function mdd = dsmh(TargetFun, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
% Dynamic Striated Metropolis-Hastings algorithm.
+% based on Waggoner/Wu/Zha (2016): "Striated Metropolis–Hastings sampler for high-dimensional models,"
+% Journal of Econometrics, 192(2): 406-420, https://doi.org/10.1016/j.jeconom.2016.02.007
%
% INPUTS
% o TargetFun [char] string specifying the name of the objective
% function (posterior kernel).
-% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
% o dataset_ data structure
% o dataset_info dataset info structure
@@ -17,23 +18,8 @@ function dsmh(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M
%
% SPECIAL REQUIREMENTS
% None.
-%
-% PARALLEL CONTEXT
-% The most computationally intensive part of this function may be executed
-% in parallel. The code suitable to be executed in
-% parallel on multi core or cluster machine (in general a 'for' cycle)
-% has been removed from this function and been placed in the posterior_sampler_core.m funtion.
-%
-% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in
-% parallel and called name_function.m and name_function_core.m and ii) a second set of functions used
-% to manage the parallel computations.
-%
-% This function was the first function to be parallelized. Later, other
-% functions have been parallelized using the same methodology.
-% Then the comments write here can be used for all the other pairs of
-% parallel functions and also for management functions.
-% Copyright © 2022-2023 Dynare Team
+% Copyright © 2022-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -50,123 +36,85 @@ function dsmh(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-opts = options_.posterior_sampler_options.dsmh;
+opts = options_.posterior_sampler_options.current_options;
+
+% Set location for the simulated particles.
+SimulationFolder = CheckPath('dsmh', M_.dname);
+%delete old stale files before creating new ones
+delete_stale_file(sprintf('%s%sparticles-*.mat', SimulationFolder,filesep))
+
+% Define prior distribution
+Prior = dprior(bayestopt_, options_.prior_trunc);
-lambda = exp(bsxfun(@minus,options_.posterior_sampler_options.dsmh.H,1:1:options_.posterior_sampler_options.dsmh.H)/(options_.posterior_sampler_options.dsmh.H-1)*log(options_.posterior_sampler_options.dsmh.lambda1));
+% Set function handle for the objective
+eval(sprintf('%s = @(x) %s(x, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, mh_bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);', 'funobj', func2str(TargetFun)));
+
+lambda = exp(bsxfun(@minus,opts.H,1:1:opts.H)/(opts.H-1)*log(opts.lambda1));
c = 0.055 ;
-MM = int64(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/10) ;
+MM = int64(opts.N*opts.G/10) ;
% Step 0: Initialization of the sampler
-[param, tlogpost_iminus1, loglik, bayestopt_] = ...
- smc_samplers_initialization(TargetFun, 'dsmh', opts.particles, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
+[param, tlogpost_iminus1, loglik] = ...
+ smc_samplers_initialization(funobj, 'dsmh', opts.particles, Prior, SimulationFolder, opts.H) ;
-ESS = zeros(options_.posterior_sampler_options.dsmh.H,1) ;
-zhat = 1 ;
+ESS = zeros(opts.H,1) ;
+zhat = 0 ;
% The DSMH starts here
-for i=2:options_.posterior_sampler_options.dsmh.H
- disp('');
- disp('Tempered iteration');
- disp(i) ;
+dprintf('#Iter. lambda ESS c Accept. rate scale resample seconds')
+for i=2:opts.H
+ resampled_particle_swarm = false;
+ t0 = tic;
% Step 1: sort the densities and compute IS weigths
- [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1,loglik,param) ;
- [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param,tlogpost_iminus1,loglik,lambda,i,zhat,ESS) ;
+ [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1, loglik,param) ;
+ [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param, tlogpost_iminus1, loglik, lambda, i, zhat, ESS) ;
+ % resample if necessary?
+ if (2*ESS(i) < opts.particles) %
+ resampled_particle_swarm = true;
+ iresample = kitagawa(weights);
+ param = param(:,iresample);
+ loglik = loglik(iresample);
+ tlogpost_i = tlogpost_i(iresample);
+ weights = ones(opts.particles, 1)/opts.particles;
+ end
% Step 2: tune c_i
- c = tune_c(TargetFun,param,tlogpost_i,lambda,i,c,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
+ [c,acpt] = tune_c(funobj, param, tlogpost_i, lambda, i, c, Omegachol, weights, mh_bounds, opts, Prior) ;
% Step 3: Metropolis step
- [param,tlogpost_iminus1,loglik] = mutation_DSMH(TargetFun,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,MM,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
+ [param,tlogpost_iminus1,loglik] = mutation_DSMH(funobj, param, tlogpost_i, tlogpost_iminus1, loglik, lambda, i, c, MM, Omegachol, weights, mh_bounds, opts, Prior) ;
+ tt = toc(t0) ;
+ if resampled_particle_swarm
+ dprintf('%3u %5.4f %9.5E %5.4f %5.4f %+5.4f %3s %5.2f', i, lambda(i), ESS(i), c, acpt, zhat, 'yes', tt)
+ else
+ dprintf('%3u %5.4f %9.5E %5.4f %5.4f %+5.4f %3s %5.2f', i, lambda(i), ESS(i), c, acpt, zhat, 'no', tt)
+ end
+ %save(sprintf('%s%sparticles-%u-%u.mat', SimulationFolder, filesep(), i, opts.H), 'param', 'tlogpost', 'loglik')
end
+tlogpost = tlogpost_iminus1 + loglik*(lambda(end)-lambda(end-1));
weights = exp(loglik*(lambda(end)-lambda(end-1)));
weights = weights/sum(weights);
-indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.particles);
-distrib_param = param(:,indx_resmpl);
-
-mean_xparam = mean(distrib_param,2);
-npar = length(xparam1);
-lb95_xparam = zeros(npar,1) ;
-ub95_xparam = zeros(npar,1) ;
-for i=1:npar
- temp = sortrows(distrib_param(i,:)') ;
- lb95_xparam(i) = temp(0.025*options_.posterior_sampler_options.dsmh.particles) ;
- ub95_xparam(i) = temp(0.975*options_.posterior_sampler_options.dsmh.particles) ;
-end
+iresample = kitagawa(weights);
+particles = param(:,iresample);
+tlogpostkernel = tlogpost(iresample);
+loglikelihood = loglik(iresample);
+save(sprintf('%s%sparameters_particles_final.mat', SimulationFolder, filesep()), 'particles', 'tlogpostkernel', 'loglikelihood')
-TeX = options_.TeX;
-
-str = sprintf(' Param. \t Lower Bound (95%%) \t Mean \t Upper Bound (95%%)');
-for l=1:npar
- name = get_the_name(l,TeX,M_,estim_params_,options_.varobs);
- str = sprintf('%s\n %s \t\t %5.4f \t\t %7.5f \t\t %5.4f', str, name, lb95_xparam(l), mean_xparam(l), ub95_xparam(l));
-end
-disp(str)
-disp('')
-
-%% Plot parameters densities
-
-if TeX
- fidTeX = fopen([M_.fname '_param_density.tex'],'w');
- fprintf(fidTeX,'%% TeX eps-loader file generated by DSMH.m (Dynare).\n');
- fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
- fprintf(fidTeX,' \n');
-end
-
-number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
-bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
-kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation.
-
-plt = 1 ;
-%for plt = 1:nbplt,
-if TeX
- NAMES = [];
- TeXNAMES = [];
-end
-hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Densities');
-for k=1:npar %min(nstar,npar-(plt-1)*nstar)
- subplot(ceil(sqrt(npar)),floor(sqrt(npar)),k)
- %kk = (plt-1)*nstar+k;
- [name,texname] = get_the_name(k,TeX,M_,estim_params_,options_.varobs);
- optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.posterior_sampler_options.dsmh.particles,bandwidth,kernel_function);
- [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
- options_.posterior_sampler_options.dsmh.particles,optimal_bandwidth,kernel_function);
- plot(density(:,1),density(:,2));
- hold on
- if TeX
- title(texname,'interpreter','latex')
- else
- title(name,'interpreter','none')
- end
- hold off
- axis tight
- drawnow
-end
-dyn_saveas(hh_fig,[ M_.fname '_param_density' int2str(plt) ],options_.nodisplay,options_.graph_format);
-if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- % TeX eps loader file
- fprintf(fidTeX,'\\begin{figure}[H]\n');
- fprintf(fidTeX,'\\centering \n');
- fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%_param_density%s}\n',min(k/floor(sqrt(npar)),1),M_.fname,int2str(plt));
- fprintf(fidTeX,'\\caption{Parameter densities based on the Dynamic Striated Metropolis-Hastings algorithm.}');
- fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
- fprintf(fidTeX,'\\end{figure}\n');
- fprintf(fidTeX,' \n');
-end
-%end
+mdd = zhat;
-function [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1,loglik,param)
+function [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1, loglik,param)
[~,indx_ord] = sortrows(tlogpost_iminus1);
tlogpost_iminus1 = tlogpost_iminus1(indx_ord);
param = param(:,indx_ord);
loglik = loglik(indx_ord);
-function [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param,tlogpost_iminus1,loglik,lambda,i,zhat,ESS)
+function [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param, tlogpost_iminus1, loglik, lambda, i, zhat, ESS)
if i==1
tlogpost_i = tlogpost_iminus1 + loglik*lambda(i);
else
tlogpost_i = tlogpost_iminus1 + loglik*(lambda(i)-lambda(i-1));
end
weights = exp(tlogpost_i-tlogpost_iminus1);
-zhat = (mean(weights))*zhat ;
+zhat = zhat + log(mean(weights)) ;
weights = weights/sum(weights);
ESS(i) = 1/sum(weights.^2);
% estimates of mean and variance
@@ -175,30 +123,30 @@ z = bsxfun(@minus,param,mu);
Omega = z*diag(weights)*z';
Omegachol = chol(Omega)';
-function c = tune_c(TargetFun,param,tlogpost_i,lambda,i,c,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
-disp('tuning c_i...');
-disp('Initial value =');
-disp(c) ;
+function [c,acpt] = tune_c(funobj, param, tlogpost_i, lambda, i, c, Omegachol, weights, mh_bounds, opts, Prior)
+% disp('tuning c_i...');
+% disp('Initial value =');
+% disp(c) ;
npar = size(param,1);
-lower_prob = (.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))^5;
-upper_prob = (.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))^(1/5);
+lower_prob = (.5*(opts.alpha0+opts.alpha1))^5;
+upper_prob = (.5*(opts.alpha0+opts.alpha1))^(1/5);
stop=0 ;
while stop==0
acpt = 0.0;
- indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.G);
+ indx_resmpl = kitagawa(weights,rand(1,1),opts.G);
param0 = param(:,indx_resmpl);
tlogpost0 = tlogpost_i(indx_resmpl);
- for j=1:options_.posterior_sampler_options.dsmh.G
- for l=1:options_.posterior_sampler_options.dsmh.K
+ for j=1:opts.G
+ for l=1:opts.K
validate = 0;
while validate == 0
candidate = param0(:,j) + sqrt(c)*Omegachol*randn(npar,1);
if all(candidate >= mh_bounds.lb) && all(candidate <= mh_bounds.ub)
- [tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,lambda(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
+ [tlogpostx,loglikx] = tempered_likelihood(funobj, candidate, lambda(i), Prior);
if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
validate = 1;
if rand(1,1)<exp(tlogpostx-tlogpost0(j)) % accept
- acpt = acpt + 1/(options_.posterior_sampler_options.dsmh.G*options_.posterior_sampler_options.dsmh.K);
+ acpt = acpt + 1/(opts.G*opts.K);
param0(:,j)= candidate;
tlogpost0(j) = tlogpostx;
end
@@ -207,29 +155,29 @@ while stop==0
end
end
end
- disp('Acceptation rate =') ;
- disp(acpt) ;
- if options_.posterior_sampler_options.dsmh.alpha0<=acpt && acpt<=options_.posterior_sampler_options.dsmh.alpha1
- disp('done!');
+ % disp('Acceptation rate =') ;
+ % disp(acpt) ;
+ if opts.alpha0<=acpt && acpt<=opts.alpha1
+ % disp('done!');
stop=1;
else
if acpt<lower_prob
c = c/5;
elseif lower_prob<=acpt && acpt<=upper_prob
- c = c*log(.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))/log(acpt);
+ c = c*log(.5*(opts.alpha0+opts.alpha1))/log(acpt);
else
c = 5*c;
end
- disp('Trying with c= ') ;
- disp(c)
+ % disp('Trying with c= ') ;
+ % disp(c)
end
end
-function [out_param,out_tlogpost_iminus1,out_loglik] = mutation_DSMH(TargetFun,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,MM,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
-indx_levels = (1:1:MM-1)*options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM;
+function [out_param,out_tlogpost_iminus1,out_loglik] = mutation_DSMH(funobj, param, tlogpost_i, tlogpost_iminus1, loglik, lambda, i, c, MM, Omegachol, weights, mh_bounds, opts, Prior)
+indx_levels = (1:1:MM-1)*opts.N*opts.G/MM;
npar = size(param,1) ;
-p = 1/(10*options_.posterior_sampler_options.dsmh.tau);
-disp('Metropolis step...');
+p = 1/(10*opts.tau);
+% disp('Metropolis step...');
% build the dynamic grid of levels
levels = [0.0;tlogpost_iminus1(indx_levels)];
% initialize the outputs
@@ -237,14 +185,14 @@ out_param = param;
out_tlogpost_iminus1 = tlogpost_i;
out_loglik = loglik;
% resample and initialize the starting groups
-indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.G);
+indx_resmpl = kitagawa(weights,rand(1,1),opts.G);
param0 = param(:,indx_resmpl);
tlogpost_iminus10 = tlogpost_iminus1(indx_resmpl);
tlogpost_i0 = tlogpost_i(indx_resmpl);
loglik0 = loglik(indx_resmpl);
% Start the Metropolis
-for l=1:options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.tau
- for j=1:options_.posterior_sampler_options.dsmh.G
+for l=1:opts.N*opts.tau
+ for j=1:opts.G
u1 = rand(1,1);
u2 = rand(1,1);
if u1<p
@@ -255,7 +203,7 @@ for l=1:options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_opt
break
end
end
- indx = floor( (k-1)*options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM+1 + u2*(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM-1) );
+ indx = floor( (k-1)*opts.N*opts.G/MM+1 + u2*(opts.N*opts.G/MM-1) );
if i==1
alp = (loglik(indx)-loglik0(j))*lambda(i);
else
@@ -272,7 +220,7 @@ for l=1:options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_opt
while validate==0
candidate = param0(:,j) + sqrt(c)*Omegachol*randn(npar,1);
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
- [tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,lambda(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
+ [tlogpostx, loglikx] = tempered_likelihood(funobj, candidate, lambda(i), Prior);
if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
validate = 1;
if u2<exp(tlogpostx-tlogpost_i0(j)) % accept
@@ -290,8 +238,8 @@ for l=1:options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_opt
end
end
end
- if mod(l,options_.posterior_sampler_options.dsmh.tau)==0
- rang = (l/options_.posterior_sampler_options.dsmh.tau-1)*options_.posterior_sampler_options.dsmh.G+1:l*options_.posterior_sampler_options.dsmh.G/options_.posterior_sampler_options.dsmh.tau;
+ if mod(l,opts.tau)==0
+ rang = (l/opts.tau-1)*opts.G+1:l*opts.G/opts.tau;
out_param(:,rang) = param0;
out_tlogpost_iminus1(rang) = tlogpost_i0;
out_loglik(rang) = loglik0;