1.5. Posterior sampling
1.5.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.5.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.5.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.5.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.5.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.5.6. Algorithm : usrsmplr Sampler
1.5.7. Tuning the proposal scale
The Metropolis-Hastings samplers (rwmh, imh) draw their proposals
from a normal or t-student distribution with covariance c*tunedCov, and
the scale c (default 1) decides how well the chain mixes. sample
takes the tuning options as a second argument:
results = sample(smplr, struct('do_tuning', true));
Option |
Default |
Meaning |
|---|---|---|
|
|
tune the scale |
|
100 |
iterations per tuning batch (both samplers) |
|
0.234 |
target acceptance rate (random walk only) |
|
0.7 |
share of the previous correction kept at each update (random walk only) |
|
0.02 |
tolerance of a batch’s acceptance rate around |
|
0.05 |
tolerance of the overall acceptance rate around |
|
4 |
tuning stops after |
Random walk (rwmh). After each batch of stepsize iterations the
scale is multiplied by a correction that moves the batch’s acceptance rate
toward alpha: a batch that accepts too little shortens the steps, one
that accepts too much lengthens them. Tuning stops once max_success+1
successive batches are within TolFun_local of alpha while the overall
rate is within TolFun_global.
Independence sampler (imh). Its proposal does not move with the chain,
so its acceptance rate measures how well the proposal covers the posterior,
and the best scale is the one that maximizes it. Shrinking c toward a
target rate, as the random walk does, makes an independence proposal
narrower than the posterior and lowers its acceptance rate further (before
this was fixed, the scale collapsed to about 1e-200 on a multimodal target).
The sampler tries the scales c0*2^k, with c0 the starting c, for
one batch of stepsize iterations each, from k=4 down to k=-4:
broad proposals first, so that the chain finds where the posterior mass lies
before the narrow ones are judged from there. A first, unscored batch at
c0*2^4 moves the chain off its start, the center of the proposal, where
it accepts nearly every proposal whatever the scale. The score of a scale is
its estimated stationary acceptance rate: the chain’s states so far, each
weighted by how long the chain stayed in it, stand for draws from the
posterior, and the proposals made at that scale for draws from the proposal.
All the scales are judged against the same states, so a narrow scale at
which the chain stalls far out scores almost nothing without dragging down
the scales tried meanwhile. While the best scale is at an edge
of the grid, the grid grows by one step (up to c0*2^20 and
c0*2^-20); a parabola in log2(c) through the best scale and its two
neighbors then refines it by at most half a step, and c stays fixed for
the rest of the run (results{1}.c reports it). alpha,
persistence_rho, the two tolerances and max_success play no role.
Two warnings can come up:
rsamplers:imhTuningFailed: no proposal was accepted at any scale of the grid. The scale stays atc0.rsamplers:imhTuningAtSearchEdge: the acceptance rate was still rising atc0*2^20(orc0*2^-20), where the search stops:tunedCovis far out of scale with the posterior.
The search takes at least ten batches (10*stepsize iterations). The draws
made during the search come from several scales: set burnin to cover it
so that the draws you keep come from the tuned proposal. The search relies
on tunedCov being about the right size (the inverse Hessian at the mode
is): a proposal hundreds of times narrower than the posterior leaves the
chain drifting at every scale of the grid, and the search ends inside it.
Checkpoints
(store_file, below) save the state of either tuner, so a resumed chain
continues the search where it stopped. The tutorial
Estimation/independence_sampler_tuning tunes the independence sampler on
the posterior of a regime-switching model.
1.5.8. Delayed acceptance: surrogate-accelerated chains
The Metropolis-Hastings samplers (rwmh, imh) accept a
da_surrogate option: a cheap approximation of the log posterior (e.g.
@(x)predict(s, x) from a surrogate fitted by emulate) that screens
every proposal before the true posterior is evaluated. Only proposals that
survive the screen pay for a true evaluation, and a second accept/reject
step keeps the chain’s stationary distribution exactly the true
posterior. See the Surrogates chapter for
the full story and a worked fs2000 example (~70% of true evaluations
skipped, ~3x wall time).
1.5.9. Checkpointing and resuming chains
Long runs should not die with a crash. The Metropolis-Hastings samplers accept three persistence options:
opts.store_file = 'mychain.mat'; % checkpoint file (MAT v7.3 = HDF5)
opts.store_every = 500; % checkpoint frequency (iterations)
opts.resume = false; % restart from the checkpoint
With store_file set, the complete chain state – draws so far, current
position, rng state, adaptive-tuning state, and the delayed-acceptance
surrogate value when active – is written every store_every iterations
(atomically: a crash mid-write cannot corrupt the previous checkpoint).
Setting resume = true restarts from the checkpoint, and because the
rng and tuning states are restored exactly, the continued chain is
bit-for-bit the chain that would have run uninterrupted. A completed
run can also be extended: resume with a larger N and the chain picks
up where it stopped. With nchain > 1, _chain<k> is appended to the
file name per chain. The tempered apt and slice samplers do not
support checkpointing yet and error early rather than silently ignoring
the option.
1.5.10. 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 rise_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.