1. Estimation

RISE estimates models by three approaches through a single entry point:

Approach	When to use it	How to call
Bayesian	Informative priors, fragile identification, or model comparison via posterior model probabilities.	`estimate(m, estim_priors = priors)`
Maximum likelihood	Data-only inference, comparison with Bayesian, or classical inference on a non-DSGE shape.	`estimate(m, estim_priors = boundsOnly)` (uniform-bounds priors)
Indirect inference	Likelihood intractable; match specific moments, IRFs, or auxiliary statistics.	`indirect_inference(m, myobjective, estim_priors = priors)`

All three share the same mode-finding optimisers, the same restriction machinery, the same Hessian / standard-error pipeline, and the same visualisation tooling. The differences are entirely in the objective and in which post-estimation outputs are meaningful.

This chapter is the modern overview of the estimation surface. The in-depth treatment is split into the sub-chapters below (restrictions, priors, posterior maximization, posterior simulation, marginal data density, processing posterior draws, indirect inference), mirroring the legacy structure.

1.8. The estimation pipeline at a glance 

A complete estimation run is four steps:

% 1. Bind data and priors to a parsed model
m = set(m, data = mydata, estim_priors = priors);

% 2. Find the posterior mode
m_est = estimate(m);

% 3. (Optional) sample the posterior
[target, x0, lb, ub, SIG] = pull_objective(m_est);
obj     = rsamplers.rwmh(target, x0, lb, ub, ...
                         struct(N = 2000, tunedCov = SIG));
results = sample(obj);

% 4. (Optional) diagnose the chain
mc = mcmc(results);
mc.summary();
mc.traceplot('phi_pi');

1.9. Dispatch on properties, not class 

A modern-toolbox distinctive worth flagging up front: the estimator chooses its likelihood path based on properties of the data and the model, never the model’s class. Concretely:

If the data has missing observations or unobserved states, the estimator runs a Kalman filter (or the regime-conditional Kalman filter, if h > 1).
If the data are complete and the model admits a closed-form posterior moment (notably BVAR with a conjugate prior), the estimator short-circuits to the closed form.
For nonlinear models the estimator runs the appropriate nonlinear filter (regime-switching unscented or particle filter).

User code that constructs a model from a new shape gets correct estimation behavior for free – there is no estimator branch keyed on class name. See Modern architecture for the broader pattern.

1.10. Priors 

A prior is a struct whose fields are parameter names and whose values are cell arrays describing the prior. RISE accepts five parametrisations for a single-parameter prior, plus a Dirichlet form for transition probabilities and a rich endogenous priors facility for behavioral constraints.

Common abbreviations used below:

start – initial value for the optimizer.
mean – mean of the prior distribution.
sd – standard deviation.
lb, ub – absolute bounds (hard truncation).
lq, uq – lower and upper quantiles.
hpp1, hpp2, hpp3, hpp4 – raw hyperparameters.
p – the struct of priors; alpha – a parameter name.

1.10.1. Uniform priors 

p.alpha = {start, lb, ub};

Internally rewritten to {start, lb, ub, 'uniform(1)'}. If you set every estimated parameter this way, you have done maximum likelihood (see Ireland, 2004) – the posterior reduces to the likelihood up to a constant.

1.10.2. Means and standard deviations 

p.alpha = {start, mean, sd, 'distribution'};

With hard truncation:

p.alpha = {start, mean, sd, 'distribution', lb, ub};

With more hyperparameters (e.g. generalized beta, t-student, truncated normal):

p.alpha = {start, mean, sd, hpp3, hpp4, 'distribution'};
p.alpha = {start, mean, sd, hpp3, hpp4, 'distribution', lb, ub};

1.10.3. Quantiles 

Means and standard deviations do not always exist (e.g. inverse gamma with certain hyperparameters; Cauchy). Quantile-based parametrisation is the universal workaround:

p.alpha = {start, lq, uq, 'distribution(probability)'};

The (probability) suffix is the mass that should fall between lq and uq – typically 0.9.

With truncation and additional hyperparameters:

p.alpha = {start, lq, uq, 'distribution(probability)', lb, ub};
p.alpha = {start, lq, uq, hpp3, hpp4, 'distribution(probability)'};
p.alpha = {start, lq, uq, hpp3, hpp4, 'distribution(probability)', lb, ub};

1.10.4. Direct hyperparameters 

The previous parametrisations are convenience wrappers around the distribution’s hyperparameters. You may set them directly:

p.alpha = {start, hpp1, hpp2, 'distribution(hpp)'};

with truncation and extra hyperparameters as for the other forms. When using this parametrisation, plot the distribution and inspect it carefully before estimating; it is the most error-prone way.

1.10.5. User-defined priors via PDF 

If you need a prior RISE does not support, supply the PDF (not the log-PDF) as a function handle plus hard bounds:

p.alpha = {start, lb, ub, function_handle};

RISE constructs the CDF, the inverse CDF, and everything else internally.

1.10.6. The t-student special case 

The t-student’s third hyperparameter is the degrees of freedom; there is no fourth, so leave it as nan or empty:

p.alpha = {start, 3, 7, 10, nan, 't(.9)'};

1.10.7. Dirichlet priors for transition matrices 

For a n x n transition matrix whose entries are named pname(i,j), declare one Dirichlet per row:

p.dirichlet_1 = {stdev(1,1), pname(1,2), mean(1,2), ...
                 pname(1,3), mean(1,3), ..., pname(1,n), mean(1,n)};

with mean(1,1) + mean(1,2) + ... + mean(1,n) = 1 (the diagonal mean is implied) and stdev(1,1) the standard deviation of the diagonal element. The _1, _2, …, _n suffix is just a label – the order is irrelevant and you only need one Dirichlet per row you actually want to put a prior on.

1.11. Endogenous (system / property / behavioral) priors 

Rather than placing priors on parameters, you may place priors on model-implied behaviors: covariances, IRFs, variance and historical decompositions, first-order solution coefficients, regime probabilities at specific dates, aggregates of any of the above. The literature calls these endogenous priors (Del Negro & Schorfheide, 2008; Christiano et al., 2011), system priors (Andrle & Benes, 2013; Andrle & Plasil, 2017, 2018), property priors (Tetlow), or behavioral priors.

The general syntax mirrors the parameter-prior cell shapes, with the parameter name replaced by an expression string:

{'expression', mean, std, 'distribution'}
{'expression', lb, ub}
{'expression', lq, uq, 'distribution(probability)'}
{'expression', hpp1, hpp2, 'distribution(hpp)'}

The expression syntax covers seven classes; regime arguments default to 1 in switching models when omitted.

1.11.1. Covariance and correlation 

'cov{var1, var2, lag}'
'cov{var1, var2, lag, regime}'
'corr{var1, var2, lag}'
'corr{var1, var2, lag, regime}'

'cov{C, Y, 1, 2}' is the covariance between C and Y at lag 1 in regime 2.

1.11.2. Impulse responses 

'irf{variable, shock, horizon}'
'irf{variable, shock, horizon, regime}'

'irf{C, EA, 3, 2}' is the response of C to shock EA at horizon 3 in regime 2.

1.11.3. Variance decomposition 

Finite horizon:

'vd{variable, shock, horizon}'
'vd{variable, shock, horizon, regime}'

Long run (set horizon = inf):

'vd{variable, shock, inf}'
'vd{variable, shock, inf, regime}'

1.11.4. Historical decomposition 

'hd{variable, shock, horizon}'

1.11.5. First-order solution coefficients 

Impact of a state variable on an endogenous variable:

'Tz{variable, stateVariable}'
'Tz{variable, stateVariable, regime}'

Examples:

'Tz{C, C{-1}}' – impact of C{-1} on C in regime 1.
'Tz{C, @sig}' – impact of the perturbation parameter.
'Tz{C, EA}' – impact of shock EA.
'Tz{C, EA{+3}, 2}' – expected impact of EA three periods ahead, in regime 2.

1.11.6. Filtered / updated / smoothed probabilities 

Regime probabilities:

'mvar{f, regime, date_range}'        % filtered
'mvar{u, regime, date_range}'        % updated
'mvar{s, regime, date_range}'        % smoothed

State combinations across chains:

'mvar{f, state1 & state2 & ... & stateN, date_range}'

'mvar{f, vol_2 & pol_1, rq(2020,3)}' – filtered probability of being simultaneously in vol_2 and pol_1 in 2020Q3.

1.11.7. Aggregations and discontinuous ranges 

sum, max, min over a date or horizon range:

'sum(irf{C, EA, 1:40, 3})'
'max(mvar{f, regime_2, rq(1970,3):rq(2024,3)})'

Discontinuous ranges:

'mvar{f, regime_2, [rq(1970,3):rq(1975,3), rq(1978,3), rq(1980,3):rq(1985,4)]}'

1.11.8. User-defined endogenous priors 

A custom function on the model:

{@(model) my_function(model), mean, std, 'distribution'}
{@(model) my_function(model), lb, ub}

For advanced cases where the prior needs filtration information, the function handle, called with one input (the model object), returns a struct:

v = struct('priors',             prior_cell, ...
           'kf_filtering_level', 0);

with kf_filtering_level in:

0 – no filtering information required;
1 – one-step-ahead (filtered);
2 – updated;
3 – smoothed.

When called with two inputs (model + filtering struct), the function returns a vector of numbers evaluating the priors.

1.12. List of supported prior distributions 

The complete catalogue, with their supports. Names with _pdf suffix (e.g. beta_pdf) are accepted as synonyms.

Name	Support	Notes
`arcsine`	\([0, 1]\)
`beta`	\([0, 1]\)	Conjugate for Bernoulli/binomial.
`dirichlet`	simplex	For transition-matrix rows; see above.
`betaprime`	\([0, \infty)\)
`cauchy`	\((-\infty, \infty)\)	No mean; quantile parametrisation only.
`exponential`	\([0, \infty)\)
`folded_normal`	\([0, \infty)\)
`gamma`	\([0, \infty)\)
`gammatau`	\([0, \infty)\)	Generalised gamma-tau.
`gumbel`	\((-\infty, \infty)\)
`half_normal`	\([0, \infty)\)
`igamma` (= `igamma2`)	\([0, \infty)\)	Use for standard deviations.
`igamma1`	\([0, \infty)\)	Inverse-gamma, shape/scale parametrisation. Use for standard deviations.
`igamma2`	\([0, \infty)\)	Use for standard deviations.
`inverse_gaussian`	\((0, \infty)\)	Use for standard deviations.
`kumaraswamy`	\((0, 1)\)	Alternative to `beta` on the unit interval.
`laplace`	\((-\infty, \infty)\)
`left_triang`	\([a, b]\)	Skewed to the left.
`levy`	\([\mu, \infty)\)
`logistic`	\((-\infty, \infty)\)
`lognormal`	\((0, \infty)\)
`normal`	\((-\infty, \infty)\)
`pareto`	\([x_m, \infty)\) with \(x_m > 0\)
`raised_cosine`	\([\mu - s, \mu + s]\) with \(s > 0\)
`sichisq`	\((0, \infty)\)	Scaled inverse \(\chi^2\). Use for standard deviations.
`t`	\((-\infty, \infty)\)	DOF is the third hyperparameter; no fourth (set `nan`).
`truncated_normal`	\([\underline{\theta}, \overline{\theta}]\)
`uniform`	\([\underline{\theta}, \overline{\theta}]\)	Implied when only the three-slot form is used.
`weibull`	\([0, \infty)\)

1.13. Visualizing priors (and posteriors)

rdist.plot plots the priors and, when given an @mcmc object, overlays the posterior densities:

rdist.plot(priors);
rdist.plot(priors, plotOpts);
rdist.plot(priors, plotOpts, mc);
rdist.plot(priors, plotOpts, mc, 'linewidth', 2);

plotOpts is a struct with the following fields:

Field	Effect
`r0c0`	Vector `[nrows, ncols]` of the figure grid. Default `[3, 3]`.
`prior_trunc`	Density-tail truncation. Default `1e-3`.
`scale_IG_trunc`	Correction of `prior_trunc` for inverse-gamma families. Default `25`.
`npoints`	Number of points for the density evaluation. Default `1000`.
`do`	Struct with logical fields `priors`, `posteriors` – whether to draw each.
`par_list`	Cell array of parameter names to restrict the plot to.

Example:

priors           = struct();
priors.beta      = {.5, 95, 5,       'beta(hpp)'};
priors.delta     = {.5, 40, 860,     'beta(hpp)'};
priors.alpha     = {.5, 360, 590,    'gamma(hpp)'};
priors.rhoz      = {.5, 80, 20,      'beta(hpp)'};
priors.sdez      = {0.1, 0.01, 0.01, 'igamma(hpp)'};

plotOpts             = struct();
plotOpts.prior_trunc = 2.6e-3;

rdist.plot(priors, plotOpts, [], struct('linewidth', 2));

1.14. Where priors live in the toolbox 

By convention, calibration and priors live together in a companion mymodel_est_params.m file with the signature function [p, priors] = mymodel_est_params(). See Model file language for why these stay out of the .rs file.

Switching parameters use their flat names with chain and regime appended (phi_pi_mon_1, mon_tp_1_2, …). Transition probabilities under a Dirichlet prior use the dirichlet_* naming pattern shown above.

1.15. Maximum likelihood (no priors)

To estimate by maximum likelihood, supply only bounds for each parameter – no mean, no distribution:

priors.alpha   = {start, lb, ub};
priors.rho     = {start, lb, ub};
priors.sigma_a = {start, lb, ub};

m_est = estimate(m, data = mydata, estim_priors = priors);

Internally RISE rewrites each entry as {start, lb, ub, 'uniform(1)'}; the prior contribution becomes a constant on the box and the posterior reduces to the (truncated) likelihood. [Ireland, 2004] is the canonical DSGE reference for this practice.

What “no prior” means in practice:

The mode is the maximum-likelihood estimator, not a posterior mode.
The Hessian-derived covariance is the standard ML covariance \(H = -\partial^2 \log L / \partial \theta^2\) at the mode.
Bayesian-only outputs are reinterpreted: the marginal data density reverts to the integrated likelihood over the box, not the Bayesian posterior probability of the model.
Posterior simulation still runs – the chain explores the likelihood-times-indicator (uniform prior). Use it for profile-likelihood-style sensitivity rather than model comparison.

Mixing approaches per parameter is allowed: some parameters can carry informative priors and others can be left with bounds-only (uniform). The dimension of the prior contribution adapts automatically.

When the likelihood is intractable – nonlinear filters that fail on your sample, or quantities of interest that depend on simulated paths – consider indirect inference instead. The optimiser suite and restriction machinery are the same.

1.16. Mode finding 

With data and priors bound, estimate runs an optimizer on the log-posterior:

m_est = estimate(m);

The default optimizer is fmincon. Pick a different one with the optimizer option:

m_est = estimate(m, optimizer = 'wcsminwel');

To pass options through to the optimizer, use the cell form {name, opts}:

opt          = optimset('fmincon');
opt.MaxIter  = 2000;
opt.TolFun   = 1e-8;
opt.Display  = 'iter';
m_est = estimate(m, optimizer = {'fmincon', opt});

1.16.1. Bundled optimizers 

Optimizer	Family	Scope	Parallel	Typical use
`fmincon`	SQP / interior-point	local	via `UseParallel`	default; smooth likelihoods.
`fminunc_bnd`	quasi-Newton (BFGS)	local	no	smooth, only-bound constraints.
`fminsearch_bnd`	Nelder-Mead	local	no	very small problems; non-smooth.
`patternsearch`	direct search (GPS / MADS)	local	via Global Opt Toolbox	non-smooth / noisy objectives.
`wcsminwel`	Sims’s csminwel + restart	local	no	robust on awkward DSGE surfaces.
`wnewrat`	Newton with numerical Hessian	local	no	mode polishing after a global pass.
`wgmhmaxlik`	MH-driven climb	local + exploration	no	pre-conditioner for multi-modal surfaces.
`bee_gate`	artificial bee colony	global (population)	yes	older global default.
`rise_lshade`, `rise_jso`, `rise_cma_es`, `rise_de`, `rise_abc`, `rise_aco`, `rise_bbo`, `rise_agsk`, `rise_mads`	`+globalopt` family	global	yes (per algorithm)	the modern global search menu.
`sa.gmhmaxlik`	simulated annealing	global	no	robust, slow; legacy.
`blockwise_optimization`	dispatcher	same as inner	same as inner	wraps any optimizer to optimize block-by-block; triggered automatically when `estim_blocks` is non-empty.

1.16.2. Common options accepted across optimizers 

Field	Meaning	Default
`MaxIter`	iteration cap	`1000`
`MaxFunEvals`	total function-eval cap	`inf`
`MaxTime`	wall-clock cap (seconds)	`inf`
`TolFun`	objective tolerance	`1e-6`
`TolX`	step tolerance	`1e-6`
`Display`	`'iter'` / `'on'` / `'off'`	`'off'`
`MaxNodes`	population size (population methods)	algorithm-specific
`ObjectiveLimit`	early-stop floor	`-inf`

Algorithm-specific hyperparameters (F, CR for rise_de; sigma0 for rise_cma_es; limit for rise_abc; …) are documented in legacy Stochastic Global Optimization.

1.16.3. Running several optimizers in turn 

A common strategy is global search to locate the basin, then a fast local method to polish:

m_est = estimate(m, ...
    optimizer = {'rise_lshade', 'wcsminwel'});

Each optimizer starts from the previous one’s best point.

1.16.4. User-defined optimizers 

optimizer may also be a function name or a function handle implementing:

[xfinal, ffinal, exitflag, H] = ...
    myOptimizer(fh, x0, lb, ub, options, varargin)

fh is the negative log-posterior to minimize; RISE guarantees lb <= x0 <= ub and a clean options struct. Outputs: xfinal (column vector in bounds), ffinal (scalar value fh(xfinal)), exitflag (MATLAB convention: 1 = converged, 0 = budget exhausted, negative = failed), H (n x n Hessian at xfinal, or [] if not available – the estimator falls back to a finite-difference pass).

Manual stopping is honoured by every +globalopt optimizer: drop a file ManualStoppingFile.txt in the working directory and the next iteration boundary returns the current best.

1.16.5. What `estimate` puts on the model 

After estimate returns:

m_est.estimation.posterior_maximization.mode
m_est.estimation.posterior_maximization.mode_stdev
m_est.estimation.posterior_maximization.hessian
m_est.estimation.posterior_maximization.vcov
m_est.estimation.posterior_maximization.log_post
m_est.estimation.posterior_maximization.log_lik
m_est.estimation.posterior_maximization.log_prior
m_est.estimation.posterior_maximization.log_marginal_data_density_laplace

Parameter names corresponding to the mode vector are fieldnames(m_est.estimation.priors) in declaration order.

1.17. Posterior simulation 

The mode is rarely the end. To explore the posterior, run an MCMC sampler. The modern API is class-based – each sampler is a class in +rsamplers constructed from a target function and bounds, then sampled with the sample method.

pull_objective lifts a sampler-ready objective out of the estimated model:

[target, x0, lb, ub, SIG] = pull_objective(m_est);

target is the log-posterior at a parameter vector. It is the value to maximize (note the sign: this differs from the legacy sampler_rwmh shortcut, which expected a value to minimize).
x0, lb, ub are the starting point and bounds.
SIG is the Hessian-derived covariance, a sensible proposal covariance.

1.17.1. Available samplers 

Class	Algorithm
`rsamplers.rwmh`	Random-walk Metropolis-Hastings. Defaults: `c=1`, `proposal='normal'`, `tunedCov` from the mode Hessian.
`rsamplers.imh`	Independent Metropolis-Hastings.
`rsamplers.apt`	Adaptive parallel tempering. Useful for multi-modal posteriors.
`rsamplers.slice`	Slice sampler. No proposal covariance to tune.
`rsamplers.usrsmplr`	User-defined sampler hook.

Options shared by every sampler (inherited from rsamplers.rsampler):

Field	Default	Role
`N`	`2000`	Draws per chain.
`nchain`	`1`	Number of chains.
`burnin`	`0`	Burn-in draws (discarded).
`thinning`	`1`	Keep every `thinning`-th draw.
`MaxTime`	`inf`	Wall-clock cap (seconds).
`MaxFunEvals`	`inf`	Total function-eval cap.

1.17.2. Constructing and running 

opts = struct( ...
    N        = 2000, ...
    nchain   = 2, ...
    burnin   = 500, ...
    c        = 0.5, ...
    tunedCov = SIG);

obj     = rsamplers.rwmh(target, x0, lb, ub, opts);
results = sample(obj);

results is a cell of length nchain. Each results{k} carries pop (struct array of draws with .x (parameter vector) and .f (log-posterior), stats, SIG, c, x0, best. Extract draws as a matrix:

draws = [results{1}.pop.x];   % nparams x N

Note

The legacy toolbox exposed convenience shortcuts sampler_rwmh, sampler_imh, sampler_apt, sampler_slice that expected a value to minimize and ran the sampler in one call. The modern path uses the class directly. The sampler_* shortcuts still work but should be migrated.

1.17.3. User-defined sampler via `rsamplers.usrsmplr`

For samplers not in the bundle (e.g. DRAM, DA-DRAM, NUTS), write a wrapper that follows the RISE output convention: a cell of nchain structs, each with the pop/stats/SIG/c fields. The legacy Posterior simulation chapter contains a worked DRAM wrapper that illustrates the pattern; it remains canonical.

1.18. Marginal data density 

For Bayesian model comparison, the mdd class wraps several estimators of \(\log p(y)\):

Laplace approximation – automatic; available without a posterior simulation.
Modified harmonic mean (Geweke).
Bridge sampling (Meng-Wong).
Mueller (importance / reciprocal-importance).
Chib-Jeliazkov (auxiliary Gibbs).

The Laplace approximation is cached on the estimated model:

log_mdd_laplace = ...
    m_est.estimation.posterior_maximization.log_marginal_data_density_laplace;

For higher-quality estimators:

md = mdd(results, m_est);          % uses the smc chain results
log_mdd_mhm = md.modified_harmonic_mean();
log_mdd_bs  = md.bridge_sampling();
log_mdd_cj  = md.chib_jeliazkov();

Each method has its own option struct; see the legacy chapter for the per-estimator option surface.

1.19. Chain diagnostics: `@mcmc`

The @mcmc class processes the output of sample and provides the standard convergence and quality diagnostics. The class itself has not changed between legacy and modern.

Constructor:

mc = mcmc(results);
mc = mcmc(results, chain_names);

Properties

Property	Content
`pnames`	Names of the estimated parameters.
`nchains`	Number of chains.
`chain_names`	Names assigned to the chains.
`npop`	Number of population (parameter) draws.
`nparams`	Number of parameters.
`draws`	MCMC parameter draws stored as a 3D array (`nparams x npop x nchains`).
`psrf`	Potential scale reduction factor (Gelman-Rubin).
`best`	Best parameter draw.
`log_post`	Log-posterior values for each draw.

Diagnostic methods

Method	Effect
`traceplot(pname)`	Per-chain trace plot of a parameter.
`densplot(pname)`	Posterior density of a parameter.
`autocorrplot(pname)`	Autocorrelation of a parameter’s chain.
`meanplot(pname)`	Running mean of a parameter’s chain.
`scatterplot(pname1, pname2)`	Joint scatter of two parameters.
`brooks_gelman()`	Brooks-Gelman PSRF (requires `nchain > 1`).
`plot_brooks_gelman()`	Plot the PSRF.
`gelman_rubin()`	Gelman-Rubin diagnostic.
`geweke(pname)`	Geweke convergence diagnostic.
`raftery_lewis(pname)`	Raftery-Lewis diagnostic.
`inefficiency_factors(pname)`	Inefficiency factors.
`summary()`	Compact tabular summary.
`posterior_parameter_statistics()`	Posterior moments + credible intervals.

Warning

mcmc is a diagnostics class constructed from existing draws. mcmc(m_est, ...) is not a sampler call and will error with “Brace indexing is not supported”. To sample, use a rsamplers.* class; to diagnose, pass the result to mcmc.

1.20. Estimation restrictions 

Linear and nonlinear restrictions on the estimated parameters are declared at the model level via the estim_linear_restrictions and estim_nonlinear_restrictions options on set / estimate. The general form is a cellstr of restrictions, one per row:

m = set(m, ...
    estim_linear_restrictions = { ...
        'phi_pi - phi_y' ; ...
        'mon_tp_1_2 - mon_tp_2_1' });

Each restriction is an expression that should equal zero. For inequalities, write them as differences with the appropriate sign convention. See the legacy chapter Estimation restrictions for the nonlinear and block-restriction extensions.

1.21. The pipeline, skeletally 

% data, parameters, priors
m     = set(m, data = db, parameters = p, estim_priors = priors);

% mode
m_est = estimate(m);

% a posterior sample
[target, x0, lb, ub, SIG] = pull_objective(m_est);
obj     = rsamplers.rwmh(target, x0, lb, ub, ...
                          struct(N = 2000, tunedCov = SIG));
results = sample(obj);

% diagnose, compute the MDD
mc          = mcmc(results);
log_mdd_mhm = mdd(results, m_est).modified_harmonic_mean();