Alternative samplers: header fixes and indentation
Johannes Pfeifer
parent
1451bf64a2
commit
13fe810453
|
|
@ -1,15 +1,11 @@
|
|||
function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||
% function DSMH_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||
% function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||
% Dynamic Striated Metropolis-Hastings algorithm.
|
||||
%
|
||||
% INPUTS
|
||||
% o TargetFun [char] string specifying the name of the objective
|
||||
% function (posterior kernel).
|
||||
% o ProposalFun [char] string specifying the name of the proposal
|
||||
% density
|
||||
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
|
||||
% o sampler_options structure
|
||||
% .invhess [double] (p*p) matrix, posterior covariance matrix (at the mode).
|
||||
% 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
|
||||
|
|
@ -37,7 +33,7 @@ function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_
|
|||
% Then the comments write here can be used for all the other pairs of
|
||||
% parallel functions and also for management functions.
|
||||
|
||||
% Copyright (C) 2006-2017 Dynare Team
|
||||
% Copyright (C) 2006-2021 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
|
|
@ -56,7 +52,7 @@ function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_
|
|||
|
||||
|
||||
lambda = exp(bsxfun(@minus,options_.dsmh.H,1:1:options_.dsmh.H)/(options_.dsmh.H-1)*log(options_.dsmh.lambda1));
|
||||
c = 0.055 ;
|
||||
c = 0.055 ;
|
||||
|
||||
% Step 0: Initialization of the sampler
|
||||
[ param, tlogpost_iminus1, loglik, npar, ~, bayestopt_] = ...
|
||||
|
|
@ -65,24 +61,24 @@ c = 0.055 ;
|
|||
ESS = zeros(options_.dsmh.H,1) ;
|
||||
zhat = 1 ;
|
||||
|
||||
% The DSMH starts here
|
||||
for i=2:options_.dsmh.H
|
||||
% The DSMH starts here
|
||||
for i=2:options_.dsmh.H
|
||||
disp('');
|
||||
disp('Tempered iteration');
|
||||
disp(i) ;
|
||||
% Step 1: sort the densities and compute IS weigths
|
||||
disp('Tempered iteration');
|
||||
disp(i) ;
|
||||
% Step 1: sort the densities and compute IS weigths
|
||||
[tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1,loglik,param) ;
|
||||
[tlogpost_i,weights,zhat,ESS,mu,Omegachol] = compute_IS_weights_and_moments(param,tlogpost_iminus1,loglik,lambda,i,zhat,ESS) ;
|
||||
% 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_);
|
||||
% Step 3: Metropolis step
|
||||
[param,tlogpost_iminus1,loglik] = mutation_DSMH(TargetFun,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
weights = exp(loglik*(lambda(end)-lambda(end-1)));
|
||||
weights = weights/sum(weights);
|
||||
indx_resmpl = smc_resampling(weights,rand(1,1),options_.dsmh.number_of_particles);
|
||||
distrib_param = param(:,indx_resmpl);
|
||||
distrib_param = param(:,indx_resmpl);
|
||||
|
||||
mean_xparam = mean(distrib_param,2);
|
||||
%mat_var_cov = bsxfun(@minus,distrib_param,mean_xparam) ;
|
||||
|
|
@ -123,134 +119,134 @@ kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform a
|
|||
|
||||
plt = 1 ;
|
||||
%for plt = 1:nbplt,
|
||||
if TeX
|
||||
NAMES = [];
|
||||
TeXNAMES = [];
|
||||
end
|
||||
hh = 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_);
|
||||
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.dsmh.number_of_particles,bandwidth,kernel_function);
|
||||
[density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
|
||||
options_.dsmh.number_of_particles,optimal_bandwidth,kernel_function);
|
||||
plot(density(:,1),density(:,2));
|
||||
hold on
|
||||
if TeX
|
||||
NAMES = [];
|
||||
TeXNAMES = [];
|
||||
end
|
||||
hh = 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_);
|
||||
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.dsmh.number_of_particles,bandwidth,kernel_function);
|
||||
[density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
|
||||
options_.dsmh.number_of_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,[ 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');
|
||||
title(texname,'interpreter','latex')
|
||||
else
|
||||
title(name,'interpreter','none')
|
||||
end
|
||||
hold off
|
||||
axis tight
|
||||
drawnow
|
||||
end
|
||||
dyn_saveas(hh,[ 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
|
||||
|
||||
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,mu,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 ;
|
||||
weights = weights/sum(weights);
|
||||
ESS(i) = 1/sum(weights.^2);
|
||||
% estimates of mean and variance
|
||||
mu = param*weights;
|
||||
z = bsxfun(@minus,param,mu);
|
||||
Omega = z*diag(weights)*z';
|
||||
Omegachol = chol(Omega)';
|
||||
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,mu,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 ;
|
||||
weights = weights/sum(weights);
|
||||
ESS(i) = 1/sum(weights.^2);
|
||||
% estimates of mean and variance
|
||||
mu = param*weights;
|
||||
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) ;
|
||||
npar = size(param,1);
|
||||
lower_prob = (.5*(options_.dsmh.alpha0+options_.dsmh.alpha1))^5;
|
||||
upper_prob = (.5*(options_.dsmh.alpha0+options_.dsmh.alpha1))^(1/5);
|
||||
stop=0 ;
|
||||
while stop==0
|
||||
acpt = 0.0;
|
||||
indx_resmpl = smc_resampling(weights,rand(1,1),options_.dsmh.G);
|
||||
param0 = param(:,indx_resmpl);
|
||||
tlogpost0 = tlogpost_i(indx_resmpl);
|
||||
for j=1:options_.dsmh.G
|
||||
for l=1:options_.dsmh.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_);
|
||||
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_.dsmh.G*options_.dsmh.K);
|
||||
param0(:,j)= candidate;
|
||||
tlogpost0(j) = tlogpostx;
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
disp('Acceptation rate =') ;
|
||||
disp(acpt) ;
|
||||
if options_.dsmh.alpha0<=acpt && acpt<=options_.dsmh.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_.dsmh.alpha0+options_.dsmh.alpha1))/log(acpt);
|
||||
else
|
||||
c = 5*c;
|
||||
disp('tuning c_i...');
|
||||
disp('Initial value =');
|
||||
disp(c) ;
|
||||
npar = size(param,1);
|
||||
lower_prob = (.5*(options_.dsmh.alpha0+options_.dsmh.alpha1))^5;
|
||||
upper_prob = (.5*(options_.dsmh.alpha0+options_.dsmh.alpha1))^(1/5);
|
||||
stop=0 ;
|
||||
while stop==0
|
||||
acpt = 0.0;
|
||||
indx_resmpl = smc_resampling(weights,rand(1,1),options_.dsmh.G);
|
||||
param0 = param(:,indx_resmpl);
|
||||
tlogpost0 = tlogpost_i(indx_resmpl);
|
||||
for j=1:options_.dsmh.G
|
||||
for l=1:options_.dsmh.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_);
|
||||
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_.dsmh.G*options_.dsmh.K);
|
||||
param0(:,j)= candidate;
|
||||
tlogpost0(j) = tlogpostx;
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
disp('Trying with c= ') ;
|
||||
disp(c)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
disp('Acceptation rate =') ;
|
||||
disp(acpt) ;
|
||||
if options_.dsmh.alpha0<=acpt && acpt<=options_.dsmh.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_.dsmh.alpha0+options_.dsmh.alpha1))/log(acpt);
|
||||
else
|
||||
c = 5*c;
|
||||
end
|
||||
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,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
|
||||
indx_levels = (1:1:options_.dsmh.M-1)*options_.dsmh.N*options_.dsmh.G/options_.dsmh.M;
|
||||
npar = size(param,1) ;
|
||||
p = 1/(10*options_.dsmh.tau);
|
||||
disp('Metropolis step...');
|
||||
% build the dynamic grid of levels
|
||||
levels = [0.0;tlogpost_iminus1(indx_levels)];
|
||||
% initialize the outputs
|
||||
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_.dsmh.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_.dsmh.N*options_.dsmh.tau
|
||||
for j=1:options_.dsmh.G
|
||||
indx_levels = (1:1:options_.dsmh.M-1)*options_.dsmh.N*options_.dsmh.G/options_.dsmh.M;
|
||||
npar = size(param,1) ;
|
||||
p = 1/(10*options_.dsmh.tau);
|
||||
disp('Metropolis step...');
|
||||
% build the dynamic grid of levels
|
||||
levels = [0.0;tlogpost_iminus1(indx_levels)];
|
||||
% initialize the outputs
|
||||
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_.dsmh.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_.dsmh.N*options_.dsmh.tau
|
||||
for j=1:options_.dsmh.G
|
||||
u1 = rand(1,1);
|
||||
u2 = rand(1,1);
|
||||
if u1<p
|
||||
|
|
@ -264,9 +260,9 @@ function [out_param,out_tlogpost_iminus1,out_loglik] = mutation_DSMH(TargetFun,p
|
|||
indx = floor( (k-1)*options_.dsmh.N*options_.dsmh.G/options_.dsmh.M+1 + u2*(options_.dsmh.N*options_.dsmh.G/options_.dsmh.M-1) );
|
||||
if i==1
|
||||
alp = (loglik(indx)-loglik0(j))*lambda(i);
|
||||
else
|
||||
else
|
||||
alp = (loglik(indx)-loglik0(j))*(lambda(i)-lambda(i-1));
|
||||
end
|
||||
end
|
||||
if u2<exp(alp)
|
||||
param0(:,j) = param(:,indx);
|
||||
tlogpost_i0(j) = tlogpost_i(indx);
|
||||
|
|
@ -274,33 +270,32 @@ function [out_param,out_tlogpost_iminus1,out_loglik] = mutation_DSMH(TargetFun,p
|
|||
tlogpost_iminus10(j) = tlogpost_iminus1(indx);
|
||||
end
|
||||
else
|
||||
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_);
|
||||
if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
|
||||
validate = 1;
|
||||
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_);
|
||||
if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
|
||||
validate = 1;
|
||||
if u2<exp(tlogpostx-tlogpost_i0(j)) % accept
|
||||
param0(:,j) = candidate;
|
||||
tlogpost_i0(j) = tlogpostx;
|
||||
param0(:,j) = candidate;
|
||||
tlogpost_i0(j) = tlogpostx;
|
||||
loglik0(j) = loglikx;
|
||||
if i==1
|
||||
tlogpost_iminus10(j) = tlogpostx-loglikx*lambda(i);
|
||||
else
|
||||
tlogpost_iminus10(j) = tlogpostx-loglikx*(lambda(i)-lambda(i-1));
|
||||
end
|
||||
if i==1
|
||||
tlogpost_iminus10(j) = tlogpostx-loglikx*lambda(i);
|
||||
else
|
||||
tlogpost_iminus10(j) = tlogpostx-loglikx*(lambda(i)-lambda(i-1));
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
if mod(l,options_.dsmh.tau)==0
|
||||
end
|
||||
if mod(l,options_.dsmh.tau)==0
|
||||
rang = (l/options_.dsmh.tau-1)*options_.dsmh.G+1:l*options_.dsmh.G/options_.dsmh.tau;
|
||||
out_param(:,rang) = param0;
|
||||
out_tlogpost_iminus1(rang) = tlogpost_i0;
|
||||
out_loglik(rang) = loglik0;
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
|
|
@ -1,15 +1,11 @@
|
|||
function Herbst_Schorfheide_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||
% function Herbst_Schorfheide_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||
% function Herbst_Schorfheide_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||
% SMC sampler from JAE 2014 .
|
||||
%
|
||||
% INPUTS
|
||||
% o TargetFun [char] string specifying the name of the objective
|
||||
% function (posterior kernel).
|
||||
% o ProposalFun [char] string specifying the name of the proposal
|
||||
% density
|
||||
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
|
||||
% o sampler_options structure
|
||||
% .invhess [double] (p*p) matrix, posterior covariance matrix (at the mode).
|
||||
% 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
|
||||
|
|
@ -37,7 +33,7 @@ function Herbst_Schorfheide_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset
|
|||
% Then the comments write here can be used for all the other pairs of
|
||||
% parallel functions and also for management functions.
|
||||
|
||||
% Copyright (C) 2006-2017 Dynare Team
|
||||
% Copyright (C) 2006-2021 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
|
|
@ -56,7 +52,7 @@ function Herbst_Schorfheide_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset
|
|||
|
||||
|
||||
% Create the tempering schedule
|
||||
phi = bsxfun(@power,(bsxfun(@minus,1:1:options_.HSsmc.nphi,1)/(options_.HSsmc.nphi-1)),options_.HSsmc.lambda) ;
|
||||
phi = bsxfun(@power,(bsxfun(@minus,1:1:options_.HSsmc.nphi,1)/(options_.HSsmc.nphi-1)),options_.HSsmc.lambda) ;
|
||||
% tuning for MH algorithms matrices
|
||||
zhat = 0 ; % normalization constant
|
||||
csim = zeros(options_.HSsmc.nphi,1) ; % scale parameter
|
||||
|
|
@ -66,7 +62,7 @@ acptsim = zeros(options_.HSsmc.nphi,1) ; % average acceptance rate
|
|||
[ param, tlogpost_i, loglik, npar, ~, bayestopt_] = ...
|
||||
SMC_samplers_initialization(TargetFun, xparam1, mh_bounds, dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,options_.HSsmc.nparticles);
|
||||
weights = ones(options_.HSsmc.nparticles,1)/options_.HSsmc.nparticles ;
|
||||
% The Herbst and Schorfheide sampler starts here
|
||||
% The Herbst and Schorfheide sampler starts here
|
||||
for i=2:options_.HSsmc.nphi
|
||||
% (a) Correction
|
||||
% incremental weights
|
||||
|
|
@ -80,26 +76,26 @@ for i=2:options_.HSsmc.nphi
|
|||
% (b) Selection
|
||||
ESSsim(i) = 1/sum(weights.^2) ;
|
||||
if (ESSsim(i) < options_.HSsmc.nparticles/2)
|
||||
indx_resmpl = smc_resampling(weights,rand(1,1),options_.HSsmc.nparticles) ;
|
||||
param = param(:,indx_resmpl) ;
|
||||
loglik = loglik(indx_resmpl) ;
|
||||
tlogpost_i = tlogpost_i(indx_resmpl) ;
|
||||
weights = ones(options_.HSsmc.nparticles,1)/options_.HSsmc.nparticles ;
|
||||
indx_resmpl = smc_resampling(weights,rand(1,1),options_.HSsmc.nparticles) ;
|
||||
param = param(:,indx_resmpl) ;
|
||||
loglik = loglik(indx_resmpl) ;
|
||||
tlogpost_i = tlogpost_i(indx_resmpl) ;
|
||||
weights = ones(options_.HSsmc.nparticles,1)/options_.HSsmc.nparticles ;
|
||||
end
|
||||
% (c) Mutation
|
||||
options_.HSsmc.c = options_.HSsmc.c*modified_logit(0.95,0.1,16.0,options_.HSsmc.acpt-options_.HSsmc.trgt) ;
|
||||
% Calculate estimates of mean and variance
|
||||
mu = param*weights ;
|
||||
z = bsxfun(@minus,param,mu) ;
|
||||
z = bsxfun(@minus,param,mu) ;
|
||||
R = z*(bsxfun(@times,z',weights)) ;
|
||||
Rchol = chol(R)' ;
|
||||
% Mutation
|
||||
% Mutation
|
||||
if options_.HSsmc.option_mutation==1
|
||||
[param,tlogpost_i,loglik,options_.HSsmc.acpt] = mutation_RW(TargetFun,param,tlogpost_i,loglik,phi,i,options_.HSsmc.c*Rchol,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_) ;
|
||||
[param,tlogpost_i,loglik,options_.HSsmc.acpt] = mutation_RW(TargetFun,param,tlogpost_i,loglik,phi,i,options_.HSsmc.c*Rchol,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_) ;
|
||||
elseif options_.HSsmc.option_mutation==2
|
||||
inv_R = inv(options_.HSsmc.c^2*R) ;
|
||||
Rdiagchol = sqrt(diag(R)) ;
|
||||
[param,tlogpost_i,loglik,options_.HSsmc.acpt] = mutation_Mixture(TargetFun,param,tlogpost_i,loglik,phi,i,options_.HSsmc.c*Rchol,options_.HSsmc.c*Rdiagchol,inv_R,mu,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_) ;
|
||||
inv_R = inv(options_.HSsmc.c^2*R) ;
|
||||
Rdiagchol = sqrt(diag(R)) ;
|
||||
[param,tlogpost_i,loglik,options_.HSsmc.acpt] = mutation_Mixture(TargetFun,param,tlogpost_i,loglik,phi,i,options_.HSsmc.c*Rchol,options_.HSsmc.c*Rdiagchol,inv_R,mu,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_) ;
|
||||
end
|
||||
acptsim(i) = options_.HSsmc.acpt ;
|
||||
csim(i) = options_.HSsmc.c ;
|
||||
|
|
@ -110,7 +106,7 @@ for i=2:options_.HSsmc.nphi
|
|||
fprintf(' %accept. = %5.4f \n', acptsim(i));
|
||||
end
|
||||
indx_resmpl = smc_resampling(weights,rand(1,1),options_.HSsmc.nparticles);
|
||||
distrib_param = param(:,indx_resmpl);
|
||||
distrib_param = param(:,indx_resmpl);
|
||||
fprintf(' Log_lik = %5.4f \n', zhat);
|
||||
|
||||
mean_xparam = mean(distrib_param,2);
|
||||
|
|
@ -151,99 +147,99 @@ 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 = 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_);
|
||||
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.HSsmc.nparticles,bandwidth,kernel_function);
|
||||
[density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
|
||||
options_.HSsmc.nparticles,optimal_bandwidth,kernel_function);
|
||||
plot(density(:,1),density(:,2));
|
||||
hold on
|
||||
if TeX
|
||||
NAMES = [];
|
||||
TeXNAMES = [];
|
||||
end
|
||||
hh = 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_);
|
||||
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.HSsmc.nparticles,bandwidth,kernel_function);
|
||||
[density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
|
||||
options_.HSsmc.nparticles,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,[ 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 Herbst/Schorfheide sampler.}');
|
||||
fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
|
||||
fprintf(fidTeX,'\\end{figure}\n');
|
||||
fprintf(fidTeX,' \n');
|
||||
title(texname,'interpreter','latex')
|
||||
else
|
||||
title(name,'interpreter','none')
|
||||
end
|
||||
hold off
|
||||
axis tight
|
||||
drawnow
|
||||
end
|
||||
dyn_saveas(hh,[ 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 Herbst/Schorfheide sampler.}');
|
||||
fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
|
||||
fprintf(fidTeX,'\\end{figure}\n');
|
||||
fprintf(fidTeX,' \n');
|
||||
end
|
||||
%end
|
||||
|
||||
function [param,tlogpost_i,loglik,acpt] = mutation_RW(TargetFun,param,tlogpost_i,loglik,phi,i,cRchol,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
|
||||
acpt = 0.0 ;
|
||||
tlogpost_i = tlogpost_i + (phi(i)-phi(i-1))*loglik ;
|
||||
for j=1:options_.HSsmc.nparticles
|
||||
validate= 0;
|
||||
while validate==0
|
||||
candidate = param(:,j) + cRchol*randn(size(param,1),1) ;
|
||||
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
|
||||
[tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,phi(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
|
||||
if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
|
||||
validate = 1;
|
||||
if rand(1,1)<exp(tlogpostx-tlogpost_i(j)) % accept
|
||||
acpt = acpt + 1 ;
|
||||
param(:,j) = candidate;
|
||||
loglik(j) = loglikx;
|
||||
tlogpost_i(j) = tlogpostx;
|
||||
end
|
||||
acpt = 0.0 ;
|
||||
tlogpost_i = tlogpost_i + (phi(i)-phi(i-1))*loglik ;
|
||||
for j=1:options_.HSsmc.nparticles
|
||||
validate= 0;
|
||||
while validate==0
|
||||
candidate = param(:,j) + cRchol*randn(size(param,1),1) ;
|
||||
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
|
||||
[tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,phi(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
|
||||
if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
|
||||
validate = 1;
|
||||
if rand(1,1)<exp(tlogpostx-tlogpost_i(j)) % accept
|
||||
acpt = acpt + 1 ;
|
||||
param(:,j) = candidate;
|
||||
loglik(j) = loglikx;
|
||||
tlogpost_i(j) = tlogpostx;
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
acpt = acpt/options_.HSsmc.nparticles;
|
||||
|
||||
end
|
||||
acpt = acpt/options_.HSsmc.nparticles;
|
||||
|
||||
function [param,tlogpost_i,loglik,acpt] = mutation_Mixture(TargetFun,param,tlogpost_i,loglik,phi,i,cRchol,cRdiagchol,invR,mu,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
|
||||
acpt = 0.0 ;
|
||||
tlogpost_i = tlogpost_i + (phi(i)-phi(i-1))*loglik ;
|
||||
for j=1:options_.HSsmc.nparticles
|
||||
qx = 0 ;
|
||||
q0 = 0 ;
|
||||
alpind = rand(1,1) ;
|
||||
validate= 0;
|
||||
while validate==0
|
||||
if alpind<options_.HSsmc.alp % RW, no need to modify qx and q0
|
||||
candidate = param(:,j) + cRchol*randn(size(param,1),1);
|
||||
elseif alpind<options_.HSsmc.alp + (1-options_.HSsmc.alp/2) % random walk with diagonal cov no need to modify qx and q0
|
||||
candidate = param(:,j) + cRdiagchol*randn(size(param,1),1);
|
||||
else % Proposal densities
|
||||
candidate = mu + cRchol*randn(size(param,1),1);
|
||||
qx = -.5*(candidate-mu)'*invR*(candidate-mu) ; % no need of the constants in the acceptation rule
|
||||
q0 = -.5*(param(:,j)-mu)'*invR*(param(:,j)-mu) ;
|
||||
end
|
||||
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
|
||||
[tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,phi(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
|
||||
if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
|
||||
validate = 1;
|
||||
if rand(1,1)<exp((tlogpostx-qx)-(tlogpost_i(j)-q0)) % accept
|
||||
acpt = acpt + 1 ;
|
||||
param(:,j) = candidate;
|
||||
loglik(j) = loglikx;
|
||||
tlogpost_i(j) = tlogpostx;
|
||||
end
|
||||
acpt = 0.0 ;
|
||||
tlogpost_i = tlogpost_i + (phi(i)-phi(i-1))*loglik ;
|
||||
for j=1:options_.HSsmc.nparticles
|
||||
qx = 0 ;
|
||||
q0 = 0 ;
|
||||
alpind = rand(1,1) ;
|
||||
validate= 0;
|
||||
while validate==0
|
||||
if alpind<options_.HSsmc.alp % RW, no need to modify qx and q0
|
||||
candidate = param(:,j) + cRchol*randn(size(param,1),1);
|
||||
elseif alpind<options_.HSsmc.alp + (1-options_.HSsmc.alp/2) % random walk with diagonal cov no need to modify qx and q0
|
||||
candidate = param(:,j) + cRdiagchol*randn(size(param,1),1);
|
||||
else % Proposal densities
|
||||
candidate = mu + cRchol*randn(size(param,1),1);
|
||||
qx = -.5*(candidate-mu)'*invR*(candidate-mu) ; % no need of the constants in the acceptation rule
|
||||
q0 = -.5*(param(:,j)-mu)'*invR*(param(:,j)-mu) ;
|
||||
end
|
||||
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
|
||||
[tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,phi(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
|
||||
if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
|
||||
validate = 1;
|
||||
if rand(1,1)<exp((tlogpostx-qx)-(tlogpost_i(j)-q0)) % accept
|
||||
acpt = acpt + 1 ;
|
||||
param(:,j) = candidate;
|
||||
loglik(j) = loglikx;
|
||||
tlogpost_i(j) = tlogpostx;
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
acpt = acpt/options_.HSsmc.nparticles;
|
||||
end
|
||||
end
|
||||
acpt = acpt/options_.HSsmc.nparticles;
|
||||
|
||||
function x = modified_logit(a,b,c,d)
|
||||
tmp = exp(c*d) ;
|
||||
x = a + b*tmp/(1 + tmp) ;
|
||||
tmp = exp(c*d) ;
|
||||
x = a + b*tmp/(1 + tmp) ;
|
||||
|
|
|
|||
Loading…
Reference in New Issue