1.5.1. The properties
pnames
mcmc/pnames is a property.
nchains
mcmc/nchains is a property.
chain_names
mcmc/chain_names is a property.
npop
mcmc/npop is a property.
nparams
mcmc/nparams is a property.
draws
mcmc/draws is a property.
psrf
mcmc/psrf is a property.
best
mcmc/best is a property.
log_post
mcmc/log_post is a property.
1.5.2. The methods
autocorrplot
Plots autocorrelations of a given parameter
hdl=autocorrplot(obj,pname) hdl=autocorrplot(obj,pname,chain_id) hdl=autocorrplot(obj,pname,chain_id,order)Args:
obj (mcmc object): mcmc object
pname (string): parameter name
- chain_id (integer | {[]}): choice of chain for which to plot the
autocorrelation. If empty, all chains are used.
order (integer | {40}): maximum order of the autocorrelation
- Returns:
:
hdl [integer]: handle to the plot
brooks_gelman
brooks_gelman computes the [Brooks and Gelman, 1998] convergence diagnostics, both the parameteric and the non-parameteric versions
INPUTS:
this : mcmc object
do_plot : [true|{false}|empty]
OUTPUTS:
myoutput [struct] : Contains
UDIAG [by 6] double
1st column: length of total sequence interval
2nd column: sum of length of within sequence intervals; used to compute mean length of within sequence intervals
3nd column: within sequence variance
4nd column: sum of within sequence variances; used to compute mean within sequence variances
5nd column: within sequence kurtosis
6nd column: sum of within sequence kurtoses; used to compute mean within sequence kurtoses
MDIAG [by 6] double
1st column: length of total sequence interval
2nd column: sum of length of within sequence intervals; used to compute mean length of within sequence intervals
3nd column: within sequence variance
4nd column: sum of within sequence variances; used to compute mean within sequence variances
5nd column: within sequence kurtosis
6nd column: sum of within sequence kurtoses; used to compute mean within sequence kurtoses
h [numeric|empty] : handle to the figures for the multivariate and univariate plots
densplot
Plots marginal density of a given parameter
hdl=densplot(obj,pname) hdl=densplot(obj,pname,chain_id) hdl=densplot(obj,pname,chain_id,N)Args:
obj (mcmc object): mcmc object
pname (string): parameter name
- chain_id (integer | {[]}): choice of chain for which to plot the
density. If empty, all chains are used.
N (integer | {250}): Number of point in the kernel density
- Returns:
:
hdl [integer]: handle to the plot
draws2struct
DRAWS2STRUCT Converts MCMC draws to a struct array of individual draws
- INPUT:
- kmcmc object with fields including
k.draws [nparams x ndraws]
k.log_post [1 x ndraws] or [ndraws x 1]
- OUTPUT:
- out1 x ndraws struct array with fields:
x : [nparams x 1] vector of parameters
f : scalar log posterior value
geweke
geweke computes Geweke convergence diagnostics
convergence=geweke(this) convergence=geweke(this,Opts) convergence=geweke(this,Opts,chain_id)Args:
this (mcmc object): mcmc object
pname (string): parameter name
chain_id (integer | {1}): choice of chain for which to do the convergence analysis
- Returns:
:
convergence [struct]: convergence analysis for each parameter
inefficiency_factors
Computes inefficiency_factors of a given parameter
Ifac=inefficiency_factors(obj) Ifac=inefficiency_factors(obj,pname) Ifac=inefficiency_factors(obj,pname,order)Args:
obj (mcmc object): mcmc object
pname (string|cell array): parameter name
- Returns:
:
Ifac [scalar|vector,struct]: inefficiency factors
mcmc
Constructor for mcmc object
mcmc_helper = mcmc(draws) mcmc_helper = mcmc(draws, pnames) mcmc_helper = mcmc(draws, pnames, subset) mcmc_helper = mcmc(draws, pnames, subset, ilinres)Args:
draws (struct): output from samplers in 3 possible forms
The draws are stored to disk in which case “draws” is a string representing the location of the draws
a structure, output of the samplers when draws are not stored to disk. In this case there is only one chain
a cell array of structures, with each structure as above. Each entry then represents one chain
a two-element cell array {draws,-logposterior} such that the first entry is a np x ndraws x nchain array of parameter values and the second entry is a 1 x ndraws x nchain array of -logPosterior values
pnames (cellstr): cell of parameter names
subset (cell array|{empty}): When not empty, subset is a 1 x 2 cell array in which the first cell contains a vector selecting the columns to retain in each chain and the second column contains the chains retained. Any or both of those cell arrays can be empty. Whenever an entry is empty, all the information available is selected. E.g. subsetting with dropping and trimming mysubs={a:b:c,[1,3,5]}. In this example, the first element selected is the one in position “a” and thereafter every “b” element is selected until we reach element in position “c”. At the same time, we select markov chains 1,3 and 5.
ilinres (function handle|{empty}): function handle that untransforms the parameters in the presence of linear restrictions.
Returns:
mcmc_helper : mcmc object
- Note:
It is the responsibility of the user to provide the names of the parameters as this routine aims to be independent from any estimation procedure or class. If the priors are set in separate structure, their names can easily be obtained via
pnames = fieldnames(priors);Alternatively, if using a RISE object, parameter names can be obtained via:
pnames = fieldnames(model.estimation.priors);
meanplot
Make a plot of cumulative means from the mcmc chain
hdl = meanplot(mcobj,pname); hdl = meanplot(mcobj,pname,chain_id);Args:
mcobj (mcmc object): mcmc object
pname (str): parameter to make the mean plot
chain_id (vector of int): id of the chain to use
- Returns:
:
hdl (handle object): handle to plot object
plot_brooks_gelman
Make a plot of cumulative means from the mcmc chain
h=plot_brooks_gelman(this,myoutput) h=plot_brooks_gelman(this,myoutput,nparamsPerFigure)Args:
this (mcmc object): mcmc object
myoutput (struct): output of brooks_gelman
nparamsPerFigure [numeric|{3}]: number of univariate parameters per figure
- Returns:
:
h (handle object): handle to all the plots (unvariates and multivariate)
posterior_parameter_statistics
posterior_parameter_statistics computes posterior parameter statistics
s = posterior_parameter_statistics(this) s = posterior_parameter_statistics(this,Opts)Args:
this [mcmc object]: mcmc object
Opts [empty|struct]: options with the following fields
percnt [vector|{[2.5,25,50,75,97.5]}] : vector used in the computation of percentiles
- -** mh_conf_sig** [scalar|{0.9}]used in the computation of Highest
Probability Density Intervals.
npoints [integer|{512}] : number of points used in the computation of the density for each parameter
kernel [char|{‘normal’}] : kernel used in the computation of the density. Choices are ‘epanechnikov’|’triangular’|’triweight’,… ‘uniform’|’cosine’
do [struct] : structure whose fields are the same as the ones for the output. Each item should be true or false, with true triggering the computation of the relevant item.
Returns:
s (struct): output with fields
posterior_mean [struct] : posterior mean
posterior_var [struct] : posterior variance
posterior_std [struct] : posterior standard deviation
posterior_percentiles [struct] : posterior posterior_percentiles as defined by the user
posterior_deciles [struct] : posterior deciles 0.1:0.1:0.9
posterior_hpd_interval [struct] : Highest probability density intervals
posterior_density [struct] : kernel density estimate
process_draws
- mcmc.process_draws is a function.
[d, fd, summary] = mcmc.process_draws(draws, subset)
psrf_plot
Makes a plot of the “posterior scale reduction factor” i.e. the [Gelman and Rubin, 1992] diagnotics from the chains.
hdl = psrf_plot(obj, pname)Args:
obj (mcmc object): mcmc object
pname (char): parameter name. N.B: One of the parameter’s name is “multivariate_” and it represents the aggregated statistics.
- start (numeric|{1}|function handle): iteration at which to start the
plot of the PSRF. If a function handle is used then it should take as input the total number of observations and return the point at which to start. e.g. @(x)round(0.5*x)
Returns:
hdl (handle object): handle to plot object
Warning:
This function requires multiple chains of MCMC samples. See nchain option of samplers.
raftery_lewis
raftery_lewis computes [Raftery and Lewis, 1992] convergence diagnostics
convergence=raftery_lewis(this) convergence=raftery_lewis(this,Opts) convergence=raftery_lewis(this,Opts,chain_id)Args:
this (mcmc object): mcmc object
pname (string): parameter name
chain_id (integer | {1}): choice of chain for which to do the convergence analysis
- Returns:
:
convergence [struct]: convergence analysis for each parameter
reload_draws
subset is a 1 x 2 cell array in which - the first cell contains the columns to retain: can be empty, defaults to all - the second column contains the chains to retain: can be empty, defaults to all
scatterplot
Make a scatter plot from the mcmc chains
hdl = scatterplot(mcobj,pname1,pname2); hdl = scatterplot(mcobj,pname1,pname2,chain_id); hdl = scatterplot(mcobj,pname1,pname2,chain_id,varargin);
- Args:
mcobj (mcmc object): mcmc object pname1 (str): x-axis parameter of scatter plot pname2 (str): y-axis parameter of scatter plot chain_id (vector of int): id of the chain to use varargin (varargin): options fed into scatter function of matlab
- Returns:
:
hdl (handle object): handle to plot object
statistics
statistics : statistics from the mcmc chains
s=statistics(obj); s=statistics(obj,opts);Args:
obj (mcmc object): mcmc object
opts (struct): structure with options of which statistics to compute. for more information see help utils.mcmc.statistics
- Returns:
:
s (struct): structure with statistics of interest
statTab (table): table representation of some of the statistics
summary
Print summary information about the mcmc draws
[myMeanStdev, MyQuantiles] = summary(mcobj); [myMeanStdev, MyQuantiles] = summary(mcobj, varargin);Args:
mcobj (mcmc object): mcmc object
varargin (options): options need to come in pairs:
‘percnt’: Quantile points (default: [2.5, 25, 50, 75, 97.5])
‘batch_size’: sampling size of for the computation of thinned standard deviation (SD(batch)), every ‘batch_size’th points are used to compute SD(batch).
- Returns:
:
myMeanStdev [cell]: mean and standard deviations
myQuantiles [cell]: quantile values
traceplot
Make a trace plot from the mcmc chains
hdl = traceplot(mcobj,pname); hdl = traceplot(mcobj,pname,chain_id); hdl = traceplot(mcobj,pname,chain_id,ma_window);Args:
mcobj (mcmc object): mcmc object
pname (str): parameter to make the trace plot
chain_id (vector of int): id of the chain to use
ma_window (int): window size if using moving average smoothing.
- Returns:
:
hdl (handle object): handle to plot object