1.4. Posterior sampling
1.4.1. The sampling function
rsamplers.rsampler.sample is the function that runs the sampling for different chains on an rsamplers.rsampler object to get a sample.
The syntax is
results = sample(obj,varargin)Args:
obj is a member of one of the following classes:
rsamplers.rwmh : Random-Walk Metropolis-Hastings
rsamplers.imh : Independent Metropolis-Hastings
rsamplers.apt : Adaptive parallel tempering
rsamplers.slice : slice sampler
rsamplers.usrsmplr : user-defined sampler
*varargin* are class-specific extra arguments
Hint
After sampling, RISE has functions for thinning, burnin, etc. through subsetting. See e.g. help on mdd (constructor for marginal data density objects).
Warning
please check the sign of the target function and make sure the algorithm used conforms with that sign. Of course, changing the sign of the target function is easy
newtarget=@(varargin)-target(varargin{:})Hint
For your convenience, some shortcuts are available so that you do not have to first create an object and then do the sampling. Moreover you do not have to create a new target if the problem is already a minimization problem. Those shortcuts are : sampler_rwmh, sampler_imh, sampler_apt, sampler_slice. They can be called using one of the following syntaxes (for the details for the inputs, see the base functions):
results = sampler_XXXX(target,x0,lb,ub)
results = sampler_XXXX(target,x0,lb,ub,opts)
1.4.2. Algorithm : Random-walk Metropolis Hastings
rsamplers.rwmh is the constructor for the Random-walk Metropolis Hastings algorithm
The syntax is
obj = rsamplers.rwmh(target,x0,lb,ub) obj = rsamplers.rwmh(target,x0,lb,ub,opts) obj = rsamplers.rwmh(target,x0,lb,ub,opts,constraints)INPUTS :
*target* : objective to MAXIMIZE
*x0* : vector of initial conditions
*lb* : lower bound
*ub* : upper bound
*opts* : options for the class. These include the properties from the superclass rsamplers.rsampler as well as the specific properties of the rwmh algorithm that can be found in the properties list below.
*constraints* (optional): - Matrix (n_constraints x 2) where each row [a, b] enforces x(a) <= x(b)
OUTPUT:
*obj* : object of class rsamplers.rwmh
OWN PROPERTIES AND DEFAULT VALUES:
*c* = 1 : cov scaling parameter
*tunedCov* : covariance matrix of the parameters
*proposal* = ‘normal’ : proposal in {‘normal’,’t-student’}
INHERITED PROPERTIES AND DEFAULT VALUES:
*nchain* = 1 : number of chains
*N* = 2000 : number of draws
*thinning* =1 : number of thinning draws
*burnin* = 0 : burnin sample
*MaxTime* = inf : Max Time
*MaxFunEvals* = inf :
See also
rsamplers.apt, rsamplers.imh, rsamplers.rsampler
1.4.3. Algorithm : Independent Metropolis Hastings
rsamplers.imh is the constructor for the Independent Metropolis-Hastings algorithm
The syntax is
obj = rsamplers.imh(target,x0,lb,ub) obj = rsamplers.imh(target,x0,lb,ub,opts) obj = rsamplers.imh(target,x0,lb,ub,opts,constraints)INPUTS :
*target* : objective to MAXIMIZE
*x0* : vector of initial conditions
*lb* : lower bound
*ub* : upper bound
*opts* : options for the class. These include the properties from the superclass rsamplers.rwmh as well as the specific properties of the imh algorithm that can be found in the properties list.
*constraints* (optional): - Matrix (n_constraints x 2) where each row [a, b] enforces x(a) <= x(b)
OUTPUT:
*obj* : object of class rsamplers.imh
OWN PROPERTIES AND DEFAULT VALUES: None
INHERITED PROPERTIES : properties from rsamplers.rwmh
See also
rsamplers.rwmh, rsamplers.apt, rsamplers.rsampler
1.4.4. Algorithm : Adaptive parallel tempering
rsamplers.apt is the constructor for the Adaptive Parallel tempering algorithm following [Miasojedow et al., 2013]
The syntax is
obj = rsamplers.apt(target,x0,lb,ub) obj = rsamplers.apt(target,x0,lb,ub,opts) obj = rsamplers.apt(target,x0,lb,ub,opts,constraints)INPUTS :
*target* : objective to MAXIMIZE
*x0* : vector of initial conditions
*lb* : lower bound
*ub* : upper bound
*opts* : options for the class. These include the properties from the superclass rsamplers.rwmh as well as the specific properties of the apt algorithm that can be found in the properties list below.
*constraints* (optional): - Matrix (n_constraints x 2) where each row [a, b] enforces x(a) <= x(b)
OUTPUT :
obj : object of class rsamplers.apt
OWN PROPERTIES AND DEFAULT VALUES:
*H* = 12 : Number of tempering stages
*alpha* = .234 : target acceptance rate
*alpha_swap* = .234 : target acceptance rate for swaps
*rwm_fixed_p* = 0 : prob of drawing from a fixed initial proposal
*rwm_exp* = 0.6 : Exponent of random-walk adaptation step size
*sw_exp_rm* = 0.6 : Exp. of temperature adaptation step size
*ram_adapt* = false : Use the robust AM adaptation
*separate_shape_adaptation* = true :separate covariance for each temperature
*fixed_temperatures* = false : fixed temperatures
INHERITED PROPERTIES : properties from rsamplers.rwmh
See also
rsamplers.rwmh, rsamplers.imh, rsamplers.rsampler
1.4.5. Algorithm : Slice Sampler
rsamplers.slice is the constructor for the slice algorithm following [Planas et al., 2015] and [Planas and Rossi, 2018]
The syntax is
obj = rsamplers.slice(target,x0,lb,ub) obj = rsamplers.slice(target,x0,lb,ub,opts) obj = rsamplers.slice(target,x0,lb,ub,opts,constraints)INPUTS:
*target* : objective to MAXIMIZE
*x0* : vector of initial conditions
*lb* : lower bound
*ub* : upper bound
*opts* : options for the class. These include the properties from the superclass rsamplers.rsampler.
*constraints* (optional): - Matrix (n_constraints x 2) where each row [a, b] enforces x(a) <= x(b)
OUTPUT:
*obj* : object of class rsamplers.rwmh
OWN PROPERTIES AND DEFAULT VALUES: None
INHERITED PROPERTIES AND DEFAULT VALUES (same as the inherited properties of rsamplers.rwmh):
*nchain* = 1 : number of chains
*N* = 2000 : number of draws
*thinning* =1 : number of thinning draws
*burnin* = 0 : burnin sample
*MaxTime* = inf : Max Time
*MaxFunEvals* = inf :
See also
rsamplers.apt, rsamplers.imh, rsamplers.rsampler, rsamplers.rwmh
1.4.6. Algorithm : usrsmplr Sampler
1.4.7. User-defined Algorithms: Another approach
Due to the fact that the RISE toolbox is written in a modular manner, a user can apply his own sampling algorithm. However, if he wants RISE to process the draws, the output of the sampling should be organized in a way that RISE understands. Here is an example in which we write a wrapping function wrapped_dramrun around the delayed-rejection adaptive Metropolis (DRAM) sampler available here .
1 function r=wrapped_dramrun(m,x0,SIG,N)
2 % wrapped_dramrun : RISE wrapper for the DRAM code
3
4 [~,lb,ub,x0_,vcov]=pull_objective(m);
5
6 [model,data,param,opts]=dramrun_inputs();
7
8 % call to the DRAM code
9 [results, chain] = dramrun(model, data, param, opts);
10
11 r=results2rise();
12
13 function [model,data,param,opts]=dramrun_inputs()
14
15 if isempty(SIG),SIG=vcov; end
16
17 if isempty(SIG),x0=x0_; end
18
19 % RISE assumes column vector, DRAM assumes row vector
20 param = struct(); param.par0 = x0.'; param.bounds=[lb,ub].';
21
22 opts=struct(); opts.qcov = SIG; opts.nsimu=N;% Number of samples
23
24 data = {};
25
26 model=struct();
27 model.ssfun = @(p)m.routines.likelihood(p.',m);
28 model.priorfun = @(p)log_prior_density(m,p.');
29
30 end
31
32 function r=results2rise()
33
34 draws=chain.';
35
36 r=struct();
37
38 for id=1:N
39 r.pop(id)=struct('x',draws(:,id),'f',nan);
40 end
41
42 r.stats.accept_ratio=results.accepted;
43 r.stats.funevals=nan; r.stats.time_elapsed=nan;
44 r.SIG=results.qcov; % covariance matrix now hidden
45 r.c=results.adascale^2;
46
47 end
48
49 end
As the code shows, there are three main articulations of the program :
Obtain from RISE and from the user the inputs that DRAM requires
Run DRAM
Transform the output of DRAM so that it can be further processed by RISE
Obtaining the DRAM inputs
Some samplers can be run directly using the log-posterior kernel, which can be obtain in RISE through the function pull_objective.
The DRAM code on the other hand, requires one function to evaluate log-prior and a separate function to evaluate the log-likelihood.
In RISE, we can get a handle to the function that evaluates the prior by using log_prior_density. We need to remember though that RISE stores the parameters vertically while DRAM stores them horizontally. Therefore, before evaluating the prior, we need to transpose the vector of parameters.
prior = @(theta) log_prior_density(model, theta.');
Where “model” is the dsge model object.
Likewise we can obtain the function computing the likelihood as follows:
likelihood = @(theta) - model.routines.likelihood(theta', model);
The other elements required by the DRAM code are easier to collect. These include : the starting vector (x0), the number of draws (n), the covariance matrix (qcov), additional inputs to the likelihood function if any (data).
Running DRAM
With the inputs collected, we can then go ahead and run the DRAM sampler.
[results, chain, s2chain] = DRAM code(model, data, param, options);
Transforming the DRAM output
If we want to process the output of DRAM through RISE, we need to transform the DRAM output in a form that RISE expects.
RISE expects each vector to be stored in a structure alongside the value of the posterior kernel. Unfortunately DRAM does not return the value of the posterior kernel. This is why we just set to “nan” all elements that cannot be obtained from the DRAM code.